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

    r2074 r2385  
    33// Corporation and others, Copyright (C) 2012, FasterCoin.  All Rights Reserved.
    44// This code is licensed under the terms of the Eclipse Public License (EPL).
    5 
    6 
    75
    86#include <cstdio>
     
    2422#include "mkl_spblas.h"
    2523#endif
    26 #if ABC_INSTRUMENT>1
     24#if ABC_INSTRUMENT > 1
    2725extern int abcPricing[20];
    2826extern int abcPricingDense[20];
     
    3735// Default Constructor
    3836//-------------------------------------------------------------------
    39 AbcMatrix::AbcMatrix ()
    40   : matrix_(NULL),
    41     model_(NULL),
    42     rowStart_(NULL),
    43     element_(NULL),
    44     column_(NULL),
    45     numberColumnBlocks_(0),
    46     numberRowBlocks_(0),
     37AbcMatrix::AbcMatrix()
     38  : matrix_(NULL)
     39  , model_(NULL)
     40  , rowStart_(NULL)
     41  , element_(NULL)
     42  , column_(NULL)
     43  , numberColumnBlocks_(0)
     44  , numberRowBlocks_(0)
     45  ,
    4746#ifdef COUNT_COPY
    48     countRealColumn_(NULL),
    49     countStartLarge_(NULL),
    50     countRow_(NULL),
    51     countElement_(NULL),
    52     smallestCount_(0),
    53     largestCount_(0),
    54 #endif
    55     startFraction_(0.0),
    56     endFraction_(1.0),
    57     savedBestDj_(0.0),
    58     originalWanted_(0),
    59     currentWanted_(0),
    60     savedBestSequence_(-1),
    61     minimumObjectsScan_(-1),
    62     minimumGoodReducedCosts_(-1)
     47  countRealColumn_(NULL)
     48  , countStartLarge_(NULL)
     49  , countRow_(NULL)
     50  , countElement_(NULL)
     51  , smallestCount_(0)
     52  , largestCount_(0)
     53  ,
     54#endif
     55  startFraction_(0.0)
     56  , endFraction_(1.0)
     57  , savedBestDj_(0.0)
     58  , originalWanted_(0)
     59  , currentWanted_(0)
     60  , savedBestSequence_(-1)
     61  , minimumObjectsScan_(-1)
     62  , minimumGoodReducedCosts_(-1)
    6363{
    6464}
     
    6767// Copy constructor
    6868//-------------------------------------------------------------------
    69 AbcMatrix::AbcMatrix (const AbcMatrix & rhs)
     69AbcMatrix::AbcMatrix(const AbcMatrix &rhs)
    7070{
    7171#ifndef COIN_SPARSE_MATRIX
     
    7474  matrix_ = new CoinPackedMatrix(*(rhs.matrix_), -0, -0);
    7575#endif
    76   model_=rhs.model_;
     76  model_ = rhs.model_;
    7777  rowStart_ = NULL;
    7878  element_ = NULL;
     
    8585#endif
    8686  numberColumnBlocks_ = rhs.numberColumnBlocks_;
    87   CoinAbcMemcpy(startColumnBlock_,rhs.startColumnBlock_,numberColumnBlocks_+1);
     87  CoinAbcMemcpy(startColumnBlock_, rhs.startColumnBlock_, numberColumnBlocks_ + 1);
    8888  numberRowBlocks_ = rhs.numberRowBlocks_;
    8989  if (numberRowBlocks_) {
    90     assert (model_);
    91     int numberRows=model_->numberRows();
     90    assert(model_);
     91    int numberRows = model_->numberRows();
    9292    int numberElements = matrix_->getNumElements();
    93     memcpy(blockStart_,rhs.blockStart_,sizeof(blockStart_));
    94     rowStart_=CoinCopyOfArray(rhs.rowStart_,numberRows*(numberRowBlocks_+2));
    95     element_=CoinCopyOfArray(rhs.element_,numberElements);
    96     column_=CoinCopyOfArray(rhs.column_,numberElements);
     93    memcpy(blockStart_, rhs.blockStart_, sizeof(blockStart_));
     94    rowStart_ = CoinCopyOfArray(rhs.rowStart_, numberRows * (numberRowBlocks_ + 2));
     95    element_ = CoinCopyOfArray(rhs.element_, numberElements);
     96    column_ = CoinCopyOfArray(rhs.column_, numberElements);
    9797#ifdef COUNT_COPY
    9898    smallestCount_ = rhs.smallestCount_;
    9999    largestCount_ = rhs.largestCount_;
    100     int numberColumns=model_->numberColumns();
    101     countRealColumn_=CoinCopyOfArray(rhs.countRealColumn_,numberColumns);
    102     memcpy(countStart_,rhs.countStart_,reinterpret_cast<char *>(&countRealColumn_)-
    103            reinterpret_cast<char *>(countStart_));
    104     int numberLarge = numberColumns-countStart_[MAX_COUNT];
    105     countStartLarge_=CoinCopyOfArray(rhs.countStartLarge_,numberLarge+1);
    106     numberElements=countStartLarge_[numberLarge];
    107     countElement_=CoinCopyOfArray(rhs.countElement_,numberElements);
    108     countRow_=CoinCopyOfArray(rhs.countRow_,numberElements);
     100    int numberColumns = model_->numberColumns();
     101    countRealColumn_ = CoinCopyOfArray(rhs.countRealColumn_, numberColumns);
     102    memcpy(countStart_, rhs.countStart_, reinterpret_cast< char * >(&countRealColumn_) - reinterpret_cast< char * >(countStart_));
     103    int numberLarge = numberColumns - countStart_[MAX_COUNT];
     104    countStartLarge_ = CoinCopyOfArray(rhs.countStartLarge_, numberLarge + 1);
     105    numberElements = countStartLarge_[numberLarge];
     106    countElement_ = CoinCopyOfArray(rhs.countElement_, numberElements);
     107    countRow_ = CoinCopyOfArray(rhs.countRow_, numberElements);
    109108#endif
    110109  }
    111110}
    112111
    113 AbcMatrix::AbcMatrix (const CoinPackedMatrix & rhs)
     112AbcMatrix::AbcMatrix(const CoinPackedMatrix &rhs)
    114113{
    115114#ifndef COIN_SPARSE_MATRIX
     
    119118#endif
    120119  matrix_->cleanMatrix();
    121   model_=NULL;
     120  model_ = NULL;
    122121  rowStart_ = NULL;
    123122  element_ = NULL;
     
    131130  largestCount_ = 0;
    132131#endif
    133   numberColumnBlocks_=1;
    134   startColumnBlock_[0]=0;
    135   startColumnBlock_[1]=0;
     132  numberColumnBlocks_ = 1;
     133  startColumnBlock_[0] = 0;
     134  startColumnBlock_[1] = 0;
    136135  numberRowBlocks_ = 0;
    137136  startFraction_ = 0;
     
    146145#ifdef ABC_SPRINT
    147146/* Subset constructor (without gaps). */
    148 AbcMatrix::AbcMatrix (const AbcMatrix & wholeMatrix,
    149                       int numberRows, const int * whichRows,
    150                       int numberColumns, const int * whichColumns)
     147AbcMatrix::AbcMatrix(const AbcMatrix &wholeMatrix,
     148  int numberRows, const int *whichRows,
     149  int numberColumns, const int *whichColumns)
    151150{
    152151#ifndef COIN_SPARSE_MATRIX
    153   matrix_ = new CoinPackedMatrix(*wholeMatrix.matrix_,numberRows,whichRows,
    154                                  numberColumns,whichColumns);
     152  matrix_ = new CoinPackedMatrix(*wholeMatrix.matrix_, numberRows, whichRows,
     153    numberColumns, whichColumns);
    155154#else
    156155  matrix_ = new CoinPackedMatrix(rhs, -0, -0);
     
    158157#endif
    159158  matrix_->cleanMatrix();
    160   model_=NULL;
     159  model_ = NULL;
    161160  rowStart_ = NULL;
    162161  element_ = NULL;
     
    170169  largestCount_ = 0;
    171170#endif
    172   numberColumnBlocks_=1;
    173   startColumnBlock_[0]=0;
    174   startColumnBlock_[1]=0;
     171  numberColumnBlocks_ = 1;
     172  startColumnBlock_[0] = 0;
     173  startColumnBlock_[1] = 0;
    175174  numberRowBlocks_ = 0;
    176175  startFraction_ = 0;
     
    187186// Destructor
    188187//-------------------------------------------------------------------
    189 AbcMatrix::~AbcMatrix ()
     188AbcMatrix::~AbcMatrix()
    190189{
    191190  delete matrix_;
    192   delete [] rowStart_;
    193   delete [] element_;
    194   delete [] column_;
     191  delete[] rowStart_;
     192  delete[] element_;
     193  delete[] column_;
    195194#ifdef COUNT_COPY
    196   delete [] countRealColumn_;
    197   delete [] countStartLarge_;
    198   delete [] countRow_;
    199   delete [] countElement_;
     195  delete[] countRealColumn_;
     196  delete[] countStartLarge_;
     197  delete[] countRow_;
     198  delete[] countElement_;
    200199#endif
    201200}
     
    205204//-------------------------------------------------------------------
    206205AbcMatrix &
    207 AbcMatrix::operator=(const AbcMatrix& rhs)
     206AbcMatrix::operator=(const AbcMatrix &rhs)
    208207{
    209208  if (this != &rhs) {
     
    214213    matrix_ = new CoinPackedMatrix(*(rhs.matrix_), -0, -0);
    215214#endif
    216     model_=rhs.model_;
    217     delete [] rowStart_;
    218     delete [] element_;
    219     delete [] column_;
     215    model_ = rhs.model_;
     216    delete[] rowStart_;
     217    delete[] element_;
     218    delete[] column_;
    220219#ifdef COUNT_COPY
    221     delete [] countRealColumn_;
    222     delete [] countStartLarge_;
    223     delete [] countRow_;
    224     delete [] countElement_;
     220    delete[] countRealColumn_;
     221    delete[] countStartLarge_;
     222    delete[] countRow_;
     223    delete[] countElement_;
    225224#endif
    226225    rowStart_ = NULL;
     
    234233#endif
    235234    numberColumnBlocks_ = rhs.numberColumnBlocks_;
    236     CoinAbcMemcpy(startColumnBlock_,rhs.startColumnBlock_,numberColumnBlocks_+1);
     235    CoinAbcMemcpy(startColumnBlock_, rhs.startColumnBlock_, numberColumnBlocks_ + 1);
    237236    numberRowBlocks_ = rhs.numberRowBlocks_;
    238237    if (numberRowBlocks_) {
    239       assert (model_);
    240       int numberRows=model_->numberRows();
     238      assert(model_);
     239      int numberRows = model_->numberRows();
    241240      int numberElements = matrix_->getNumElements();
    242       memcpy(blockStart_,rhs.blockStart_,sizeof(blockStart_));
    243       rowStart_=CoinCopyOfArray(rhs.rowStart_,numberRows*(numberRowBlocks_+2));
    244       element_=CoinCopyOfArray(rhs.element_,numberElements);
    245       column_=CoinCopyOfArray(rhs.column_,numberElements);
     241      memcpy(blockStart_, rhs.blockStart_, sizeof(blockStart_));
     242      rowStart_ = CoinCopyOfArray(rhs.rowStart_, numberRows * (numberRowBlocks_ + 2));
     243      element_ = CoinCopyOfArray(rhs.element_, numberElements);
     244      column_ = CoinCopyOfArray(rhs.column_, numberElements);
    246245#ifdef COUNT_COPY
    247246      smallestCount_ = rhs.smallestCount_;
    248247      largestCount_ = rhs.largestCount_;
    249       int numberColumns=model_->numberColumns();
    250       countRealColumn_=CoinCopyOfArray(rhs.countRealColumn_,numberColumns);
    251       memcpy(countStart_,rhs.countStart_,reinterpret_cast<char *>(&countRealColumn_)-
    252              reinterpret_cast<char *>(countStart_));
    253       int numberLarge = numberColumns-countStart_[MAX_COUNT];
    254       countStartLarge_=CoinCopyOfArray(rhs.countStartLarge_,numberLarge+1);
    255       numberElements=countStartLarge_[numberLarge];
    256       countElement_=CoinCopyOfArray(rhs.countElement_,numberElements);
    257       countRow_=CoinCopyOfArray(rhs.countRow_,numberElements);
     248      int numberColumns = model_->numberColumns();
     249      countRealColumn_ = CoinCopyOfArray(rhs.countRealColumn_, numberColumns);
     250      memcpy(countStart_, rhs.countStart_, reinterpret_cast< char * >(&countRealColumn_) - reinterpret_cast< char * >(countStart_));
     251      int numberLarge = numberColumns - countStart_[MAX_COUNT];
     252      countStartLarge_ = CoinCopyOfArray(rhs.countStartLarge_, numberLarge + 1);
     253      numberElements = countStartLarge_[numberLarge];
     254      countElement_ = CoinCopyOfArray(rhs.countElement_, numberElements);
     255      countRow_ = CoinCopyOfArray(rhs.countRow_, numberElements);
    258256#endif
    259257    }
     
    270268}
    271269// Sets model
    272 void
    273 AbcMatrix::setModel(AbcSimplex * model)
    274 {
    275   model_=model;
    276   int numberColumns=model_->numberColumns();
    277   bool needExtension=numberColumns>matrix_->getNumCols();
     270void AbcMatrix::setModel(AbcSimplex *model)
     271{
     272  model_ = model;
     273  int numberColumns = model_->numberColumns();
     274  bool needExtension = numberColumns > matrix_->getNumCols();
    278275  if (needExtension) {
    279276    CoinBigIndex lastElement = matrix_->getNumElements();
    280     matrix_->reserve(numberColumns,lastElement,true);
    281     CoinBigIndex * columnStart = matrix_->getMutableVectorStarts();
    282     for (int i=numberColumns;i>=0;i--) {
    283       if (columnStart[i]==0)
    284         columnStart[i]=lastElement;
     277    matrix_->reserve(numberColumns, lastElement, true);
     278    CoinBigIndex *columnStart = matrix_->getMutableVectorStarts();
     279    for (int i = numberColumns; i >= 0; i--) {
     280      if (columnStart[i] == 0)
     281        columnStart[i] = lastElement;
    285282      else
    286         break;
    287     }
    288     assert (lastElement==columnStart[numberColumns]);
     283        break;
     284    }
     285    assert(lastElement == columnStart[numberColumns]);
    289286  }
    290287}
     
    293290AbcMatrix::reverseOrderedCopy() const
    294291{
    295   CoinPackedMatrix * matrix = new CoinPackedMatrix();
     292  CoinPackedMatrix *matrix = new CoinPackedMatrix();
    296293  matrix->setExtraGap(0.0);
    297294  matrix->setExtraMajor(0.0);
     
    301298/// returns number of elements in column part of basis,
    302299CoinBigIndex
    303 AbcMatrix::countBasis( const int * whichColumn,
    304                        int & numberColumnBasic)
    305 {
    306   const int *  COIN_RESTRICT columnLength = matrix_->getVectorLengths();
     300AbcMatrix::countBasis(const int *whichColumn,
     301  int &numberColumnBasic)
     302{
     303  const int *COIN_RESTRICT columnLength = matrix_->getVectorLengths();
    307304  CoinBigIndex numberElements = 0;
    308   int numberRows=model_->numberRows();
     305  int numberRows = model_->numberRows();
    309306  // just count - can be over so ignore zero problem
    310307  for (int i = 0; i < numberColumnBasic; i++) {
    311     int iColumn = whichColumn[i]-numberRows;
     308    int iColumn = whichColumn[i] - numberRows;
    312309    numberElements += columnLength[iColumn];
    313310  }
    314311  return numberElements;
    315312}
    316 void
    317 AbcMatrix::fillBasis(const int * COIN_RESTRICT whichColumn,
    318                      int & numberColumnBasic,
    319                      int * COIN_RESTRICT indexRowU,
    320                      int * COIN_RESTRICT start,
    321                      int * COIN_RESTRICT rowCount,
    322                      int * COIN_RESTRICT columnCount,
    323                      CoinFactorizationDouble * COIN_RESTRICT elementU)
    324 {
    325   const int * COIN_RESTRICT columnLength = matrix_->getVectorLengths();
     313void AbcMatrix::fillBasis(const int *COIN_RESTRICT whichColumn,
     314  int &numberColumnBasic,
     315  int *COIN_RESTRICT indexRowU,
     316  int *COIN_RESTRICT start,
     317  int *COIN_RESTRICT rowCount,
     318  int *COIN_RESTRICT columnCount,
     319  CoinFactorizationDouble *COIN_RESTRICT elementU)
     320{
     321  const int *COIN_RESTRICT columnLength = matrix_->getVectorLengths();
    326322  CoinBigIndex numberElements = start[0];
    327323  // fill
    328   const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
    329   const int * COIN_RESTRICT row = matrix_->getIndices();
    330   const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
    331   int numberRows=model_->numberRows();
     324  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     325  const int *COIN_RESTRICT row = matrix_->getIndices();
     326  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
     327  int numberRows = model_->numberRows();
    332328  for (int i = 0; i < numberColumnBasic; i++) {
    333     int iColumn = whichColumn[i]-numberRows;
     329    int iColumn = whichColumn[i] - numberRows;
    334330    int length = columnLength[iColumn];
    335331    CoinBigIndex startThis = columnStart[iColumn];
     
    340336      indexRowU[numberElements] = iRow;
    341337      rowCount[iRow]++;
    342       assert (elementByColumn[j]);
     338      assert(elementByColumn[j]);
    343339      elementU[numberElements++] = elementByColumn[j];
    344340    }
    345     start[i+1] = numberElements;
     341    start[i + 1] = numberElements;
    346342  }
    347343}
    348344#ifdef ABC_LONG_FACTORIZATION
    349 void
    350 AbcMatrix::fillBasis(const int * COIN_RESTRICT whichColumn,
    351                      int & numberColumnBasic,
    352                      int * COIN_RESTRICT indexRowU,
    353                      int * COIN_RESTRICT start,
    354                      int * COIN_RESTRICT rowCount,
    355                      int * COIN_RESTRICT columnCount,
    356                      long double * COIN_RESTRICT elementU)
    357 {
    358   const int * COIN_RESTRICT columnLength = matrix_->getVectorLengths();
     345void AbcMatrix::fillBasis(const int *COIN_RESTRICT whichColumn,
     346  int &numberColumnBasic,
     347  int *COIN_RESTRICT indexRowU,
     348  int *COIN_RESTRICT start,
     349  int *COIN_RESTRICT rowCount,
     350  int *COIN_RESTRICT columnCount,
     351  long double *COIN_RESTRICT elementU)
     352{
     353  const int *COIN_RESTRICT columnLength = matrix_->getVectorLengths();
    359354  CoinBigIndex numberElements = start[0];
    360355  // fill
    361   const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
    362   const int * COIN_RESTRICT row = matrix_->getIndices();
    363   const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
    364   int numberRows=model_->numberRows();
     356  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     357  const int *COIN_RESTRICT row = matrix_->getIndices();
     358  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
     359  int numberRows = model_->numberRows();
    365360  for (int i = 0; i < numberColumnBasic; i++) {
    366     int iColumn = whichColumn[i]-numberRows;
     361    int iColumn = whichColumn[i] - numberRows;
    367362    int length = columnLength[iColumn];
    368363    CoinBigIndex startThis = columnStart[iColumn];
     
    373368      indexRowU[numberElements] = iRow;
    374369      rowCount[iRow]++;
    375       assert (elementByColumn[j]);
     370      assert(elementByColumn[j]);
    376371      elementU[numberElements++] = elementByColumn[j];
    377372    }
    378     start[i+1] = numberElements;
    379   }
    380 }
    381 #endif 
     373    start[i + 1] = numberElements;
     374  }
     375}
     376#endif
    382377#if 0
    383378/// Move largest in column to beginning
     
    416411#endif
    417412// Creates row copy
    418 void
    419 AbcMatrix::createRowCopy()
     413void AbcMatrix::createRowCopy()
    420414{
    421415#if ABC_PARALLEL
    422   if (model_->parallelMode()==0)
    423 #endif
    424     numberRowBlocks_=1;
     416  if (model_->parallelMode() == 0)
     417#endif
     418    numberRowBlocks_ = 1;
    425419#if ABC_PARALLEL
    426420  else
    427     numberRowBlocks_=CoinMin(NUMBER_ROW_BLOCKS,model_->numberCpus());
    428 #endif
    429   int maximumRows=model_->maximumAbcNumberRows();
    430   int numberRows=model_->numberRows();
    431   int numberColumns=model_->numberColumns();
     421    numberRowBlocks_ = CoinMin(NUMBER_ROW_BLOCKS, model_->numberCpus());
     422#endif
     423  int maximumRows = model_->maximumAbcNumberRows();
     424  int numberRows = model_->numberRows();
     425  int numberColumns = model_->numberColumns();
    432426  int numberElements = matrix_->getNumElements();
    433   assert (!rowStart_);
    434   char * whichBlock_=new char [numberColumns];
    435   rowStart_=new CoinBigIndex[numberRows*(numberRowBlocks_+2)];
    436   element_ = new double [numberElements];
    437   column_ = new int [numberElements];
    438   const CoinBigIndex *  COIN_RESTRICT columnStart = matrix_->getVectorStarts();
    439   memset(blockStart_,0,sizeof(blockStart_));
     427  assert(!rowStart_);
     428  char *whichBlock_ = new char[numberColumns];
     429  rowStart_ = new CoinBigIndex[numberRows * (numberRowBlocks_ + 2)];
     430  element_ = new double[numberElements];
     431  column_ = new int[numberElements];
     432  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     433  memset(blockStart_, 0, sizeof(blockStart_));
    440434  int ecount[10];
    441   assert (numberRowBlocks_<16);
    442   CoinAbcMemset0(ecount,10);
     435  assert(numberRowBlocks_ < 16);
     436  CoinAbcMemset0(ecount, 10);
    443437  // allocate to blocks (put a bit less in first as will be dealing with slacks) LATER
    444   CoinBigIndex start=0;
    445   int block=0;
    446   CoinBigIndex work=(2*numberColumns+matrix_->getNumElements()+numberRowBlocks_-1)/numberRowBlocks_;
    447   CoinBigIndex thisWork=work;
    448   for (int iColumn=0;iColumn<numberColumns;iColumn++) {
    449     CoinBigIndex end=columnStart[iColumn+1];
    450     assert (block>=0&&block<numberRowBlocks_);
    451     whichBlock_[iColumn]=static_cast<char>(block);
    452     thisWork -= 2+end-start;
    453     ecount[block]+=end-start;
    454     start=end;
     438  CoinBigIndex start = 0;
     439  int block = 0;
     440  CoinBigIndex work = (2 * numberColumns + matrix_->getNumElements() + numberRowBlocks_ - 1) / numberRowBlocks_;
     441  CoinBigIndex thisWork = work;
     442  for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     443    CoinBigIndex end = columnStart[iColumn + 1];
     444    assert(block >= 0 && block < numberRowBlocks_);
     445    whichBlock_[iColumn] = static_cast< char >(block);
     446    thisWork -= 2 + end - start;
     447    ecount[block] += end - start;
     448    start = end;
    455449    blockStart_[block]++;
    456     if (thisWork<=0) {
     450    if (thisWork <= 0) {
    457451      block++;
    458       thisWork=work;
     452      thisWork = work;
    459453    }
    460454  }
     
    465459  printf("\n");
    466460#endif
    467   start=0;
    468   for (int i=0;i<numberRowBlocks_;i++) {
    469     int next=start+blockStart_[i];
    470     blockStart_[i]=start;
    471     start=next;
    472   }
    473   blockStart_[numberRowBlocks_]=start;
    474   assert (start==numberColumns);
    475   const int *  COIN_RESTRICT row = matrix_->getIndices();
    476   const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
     461  start = 0;
     462  for (int i = 0; i < numberRowBlocks_; i++) {
     463    int next = start + blockStart_[i];
     464    blockStart_[i] = start;
     465    start = next;
     466  }
     467  blockStart_[numberRowBlocks_] = start;
     468  assert(start == numberColumns);
     469  const int *COIN_RESTRICT row = matrix_->getIndices();
     470  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
    477471  // counts
    478   CoinAbcMemset0(rowStart_,numberRows*(numberRowBlocks_+2));
    479   int * COIN_RESTRICT last = rowStart_+numberRows*(numberRowBlocks_+1);
    480   for (int iColumn=0;iColumn<numberColumns;iColumn++) {
    481     int block=whichBlock_[iColumn];
     472  CoinAbcMemset0(rowStart_, numberRows * (numberRowBlocks_ + 2));
     473  int *COIN_RESTRICT last = rowStart_ + numberRows * (numberRowBlocks_ + 1);
     474  for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     475    int block = whichBlock_[iColumn];
    482476    blockStart_[block]++;
    483     int base=(block+1)*numberRows;
    484     for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn+1];j++) {
    485       int iRow=row[j];
    486       rowStart_[base+iRow]++;
     477    int base = (block + 1) * numberRows;
     478    for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn + 1]; j++) {
     479      int iRow = row[j];
     480      rowStart_[base + iRow]++;
    487481      last[iRow]++;
    488482    }
    489483  }
    490484  // back to real starts
    491   for (int iBlock=numberRowBlocks_-1;iBlock>=0;iBlock--)
    492     blockStart_[iBlock+1]=blockStart_[iBlock];
    493   blockStart_[0]=0;
    494   CoinBigIndex put=0;
    495   for (int iRow=1;iRow<numberRows;iRow++) {
    496     int n=last[iRow-1];
    497     last[iRow-1]=put;
     485  for (int iBlock = numberRowBlocks_ - 1; iBlock >= 0; iBlock--)
     486    blockStart_[iBlock + 1] = blockStart_[iBlock];
     487  blockStart_[0] = 0;
     488  CoinBigIndex put = 0;
     489  for (int iRow = 1; iRow < numberRows; iRow++) {
     490    int n = last[iRow - 1];
     491    last[iRow - 1] = put;
    498492    put += n;
    499     rowStart_[iRow]=put;
    500   }
    501   int n=last[numberRows-1];
    502   last[numberRows-1]=put;
     493    rowStart_[iRow] = put;
     494  }
     495  int n = last[numberRows - 1];
     496  last[numberRows - 1] = put;
    503497  put += n;
    504   assert (put==numberElements);
     498  assert(put == numberElements);
    505499  //last[numberRows-1]=put;
    506500  // starts
    507   int base=0;
    508   for (int iBlock=0;iBlock<numberRowBlocks_;iBlock++) {
    509     for (int iRow=0;iRow<numberRows;iRow++) {
    510       rowStart_[base+numberRows+iRow]+=rowStart_[base+iRow];
    511     }
    512     base+= numberRows;
    513   }
    514   for (int iColumn=0;iColumn<numberColumns;iColumn++) {
    515     int block=whichBlock_[iColumn];
    516     int where=blockStart_[block];
     501  int base = 0;
     502  for (int iBlock = 0; iBlock < numberRowBlocks_; iBlock++) {
     503    for (int iRow = 0; iRow < numberRows; iRow++) {
     504      rowStart_[base + numberRows + iRow] += rowStart_[base + iRow];
     505    }
     506    base += numberRows;
     507  }
     508  for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     509    int block = whichBlock_[iColumn];
     510    int where = blockStart_[block];
    517511    blockStart_[block]++;
    518     int base=block*numberRows;
    519     for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn+1];j++) {
    520       int iRow=row[j];
    521       CoinBigIndex put=rowStart_[base+iRow];
    522       rowStart_[base+iRow]++;
    523       column_[put]=where+maximumRows;
    524       element_[put]=elementByColumn[j];
     512    int base = block * numberRows;
     513    for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn + 1]; j++) {
     514      int iRow = row[j];
     515      CoinBigIndex put = rowStart_[base + iRow];
     516      rowStart_[base + iRow]++;
     517      column_[put] = where + maximumRows;
     518      element_[put] = elementByColumn[j];
    525519    }
    526520  }
    527521  // back to real starts etc
    528   base=numberRows*numberRowBlocks_;
    529   for (int iBlock=numberRowBlocks_-1;iBlock>=0;iBlock--) {
    530     blockStart_[iBlock+1]=blockStart_[iBlock]+maximumRows;
    531     CoinAbcMemcpy(rowStart_+base,rowStart_+base-numberRows,numberRows);
     522  base = numberRows * numberRowBlocks_;
     523  for (int iBlock = numberRowBlocks_ - 1; iBlock >= 0; iBlock--) {
     524    blockStart_[iBlock + 1] = blockStart_[iBlock] + maximumRows;
     525    CoinAbcMemcpy(rowStart_ + base, rowStart_ + base - numberRows, numberRows);
    532526    base -= numberRows;
    533527  }
    534   blockStart_[0]=0;//maximumRows;
    535   delete [] whichBlock_;
    536   // and move 
    537   CoinAbcMemcpy(rowStart_,last,numberRows);
     528  blockStart_[0] = 0; //maximumRows;
     529  delete[] whichBlock_;
     530  // and move
     531  CoinAbcMemcpy(rowStart_, last, numberRows);
    538532  // All in useful
    539   CoinAbcMemcpy(rowStart_+(numberRowBlocks_+1)*numberRows,
    540                 rowStart_+(numberRowBlocks_)*numberRows,numberRows);
     533  CoinAbcMemcpy(rowStart_ + (numberRowBlocks_ + 1) * numberRows,
     534    rowStart_ + (numberRowBlocks_)*numberRows, numberRows);
    541535#ifdef COUNT_COPY
    542536  // now blocked by element count
    543   countRealColumn_=new int [numberColumns];
    544   int counts [2*MAX_COUNT];
    545   memset(counts,0,sizeof(counts));
     537  countRealColumn_ = new int[numberColumns];
     538  int counts[2 * MAX_COUNT];
     539  memset(counts, 0, sizeof(counts));
    546540  //memset(countFirst_,0,sizeof(countFirst_));
    547   int numberLarge=0;
    548   for (int iColumn=0;iColumn<numberColumns;iColumn++) {
    549     int n=columnStart[iColumn+1]-columnStart[iColumn];
    550     if (n<MAX_COUNT)
     541  int numberLarge = 0;
     542  for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     543    int n = columnStart[iColumn + 1] - columnStart[iColumn];
     544    if (n < MAX_COUNT)
    551545      counts[n]++;
    552546    else
    553547      numberLarge++;
    554548  }
    555   CoinBigIndex numberExtra=3; // for alignment
     549  CoinBigIndex numberExtra = 3; // for alignment
    556550#define SMALL_COUNT 1
    557   for (int i=0;i<MAX_COUNT;i++) {
    558     int n=counts[i];
    559     if (n>=SMALL_COUNT) {
     551  for (int i = 0; i < MAX_COUNT; i++) {
     552    int n = counts[i];
     553    if (n >= SMALL_COUNT) {
    560554      n &= 3;
    561       int extra=(4-n)&3;
    562       numberExtra+= i*extra;
     555      int extra = (4 - n) & 3;
     556      numberExtra += i * extra;
    563557    } else {
    564558      // treat as large
    565       numberLarge+=n;
    566     }
    567   }
    568   countElement_= new double [numberElements+numberExtra];
    569   countRow_=new int [numberElements+numberExtra];
    570   countStartLarge_ = new CoinBigIndex [numberLarge+1];
    571   countStartLarge_[numberLarge]=numberElements+numberExtra;
     559      numberLarge += n;
     560    }
     561  }
     562  countElement_ = new double[numberElements + numberExtra];
     563  countRow_ = new int[numberElements + numberExtra];
     564  countStartLarge_ = new CoinBigIndex[numberLarge + 1];
     565  countStartLarge_[numberLarge] = numberElements + numberExtra;
    572566  //return;
    573   CoinInt64 xx = reinterpret_cast<CoinInt64>(countElement_);
    574   int iBottom = static_cast<int>(xx & 31);
    575   int offset = iBottom>>3;
    576   CoinBigIndex firstElementLarge=0;
     567  CoinInt64 xx = reinterpret_cast< CoinInt64 >(countElement_);
     568  int iBottom = static_cast< int >(xx & 31);
     569  int offset = iBottom >> 3;
     570  CoinBigIndex firstElementLarge = 0;
    577571  if (offset)
    578     firstElementLarge=4-offset;
     572    firstElementLarge = 4 - offset;
    579573  //countStart_[0]=firstElementLarge;
    580   int positionLarge=0;
    581   smallestCount_=0;
    582   largestCount_=0;
    583   for (int i=0;i<MAX_COUNT;i++) {
    584     int n=counts[i];
    585     countFirst_[i]=positionLarge;
    586     countStart_[i]=firstElementLarge;
    587     if (n>=SMALL_COUNT) {
    588       counts[i+MAX_COUNT]=1;
    589       if (smallestCount_==0)
    590         smallestCount_=i;
    591       largestCount_=i;
    592       positionLarge+=n;
    593       firstElementLarge+=n*i;
     574  int positionLarge = 0;
     575  smallestCount_ = 0;
     576  largestCount_ = 0;
     577  for (int i = 0; i < MAX_COUNT; i++) {
     578    int n = counts[i];
     579    countFirst_[i] = positionLarge;
     580    countStart_[i] = firstElementLarge;
     581    if (n >= SMALL_COUNT) {
     582      counts[i + MAX_COUNT] = 1;
     583      if (smallestCount_ == 0)
     584        smallestCount_ = i;
     585      largestCount_ = i;
     586      positionLarge += n;
     587      firstElementLarge += n * i;
    594588      n &= 3;
    595       int extra=(4-n)&3;
    596       firstElementLarge+= i*extra;
    597     }
    598     counts[i]=0;
     589      int extra = (4 - n) & 3;
     590      firstElementLarge += i * extra;
     591    }
     592    counts[i] = 0;
    599593  }
    600594  largestCount_++;
    601   countFirst_[MAX_COUNT]=positionLarge;
    602   countStart_[MAX_COUNT]=firstElementLarge;
    603   numberLarge=0;
    604   for (int iColumn=0;iColumn<numberColumns;iColumn++) {
    605     CoinBigIndex start=columnStart[iColumn];
    606     int n=columnStart[iColumn+1]-start;
     595  countFirst_[MAX_COUNT] = positionLarge;
     596  countStart_[MAX_COUNT] = firstElementLarge;
     597  numberLarge = 0;
     598  for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     599    CoinBigIndex start = columnStart[iColumn];
     600    int n = columnStart[iColumn + 1] - start;
    607601    CoinBigIndex put;
    608602    int position;
    609     if (n<MAX_COUNT&&counts[MAX_COUNT+n]!=0) {
    610       int iWhich=counts[n];
    611       position = countFirst_[n]+iWhich;
    612       int iBlock4=iWhich&(~3);
     603    if (n < MAX_COUNT && counts[MAX_COUNT + n] != 0) {
     604      int iWhich = counts[n];
     605      position = countFirst_[n] + iWhich;
     606      int iBlock4 = iWhich & (~3);
    613607      iWhich &= 3;
    614       put = countStart_[n]+iBlock4*n;
     608      put = countStart_[n] + iBlock4 * n;
    615609      put += iWhich;
    616610      counts[n]++;
    617       for (int i=0;i<n;i++) {
    618         countRow_[put]=row[start+i];
    619         countElement_[put]=elementByColumn[start+i];
    620         put+=4;
     611      for (int i = 0; i < n; i++) {
     612        countRow_[put] = row[start + i];
     613        countElement_[put] = elementByColumn[start + i];
     614        put += 4;
    621615      }
    622616    } else {
    623       position = positionLarge+numberLarge;
    624       put=firstElementLarge;
    625       countStartLarge_[numberLarge]=put;
    626       firstElementLarge+=n;
     617      position = positionLarge + numberLarge;
     618      put = firstElementLarge;
     619      countStartLarge_[numberLarge] = put;
     620      firstElementLarge += n;
    627621      numberLarge++;
    628       CoinAbcMemcpy(countRow_+put,row+start,n);
    629       CoinAbcMemcpy(countElement_+put,elementByColumn+start,n);
    630     }
    631     countRealColumn_[position]=iColumn+maximumRows;
    632   }
    633   countStartLarge_[numberLarge]=firstElementLarge;
    634   assert (firstElementLarge<=numberElements+numberExtra);
     622      CoinAbcMemcpy(countRow_ + put, row + start, n);
     623      CoinAbcMemcpy(countElement_ + put, elementByColumn + start, n);
     624    }
     625    countRealColumn_[position] = iColumn + maximumRows;
     626  }
     627  countStartLarge_[numberLarge] = firstElementLarge;
     628  assert(firstElementLarge <= numberElements + numberExtra);
    635629#endif
    636630}
    637631// Make all useful
    638 void
    639 AbcMatrix::makeAllUseful(CoinIndexedVector & /*spare*/)
    640 {
    641   int numberRows=model_->numberRows();
    642   CoinBigIndex * COIN_RESTRICT rowEnd = rowStart_+numberRows;
    643   const CoinBigIndex * COIN_RESTRICT rowReallyEnd = rowStart_+2*numberRows;
    644   for (int iRow=0;iRow<numberRows;iRow++) {
    645     rowEnd[iRow]=rowReallyEnd[iRow];
     632void AbcMatrix::makeAllUseful(CoinIndexedVector & /*spare*/)
     633{
     634  int numberRows = model_->numberRows();
     635  CoinBigIndex *COIN_RESTRICT rowEnd = rowStart_ + numberRows;
     636  const CoinBigIndex *COIN_RESTRICT rowReallyEnd = rowStart_ + 2 * numberRows;
     637  for (int iRow = 0; iRow < numberRows; iRow++) {
     638    rowEnd[iRow] = rowReallyEnd[iRow];
    646639  }
    647640}
    648641#define SQRT_ARRAY
    649642// Creates scales for column copy (rowCopy in model may be modified)
    650 void
    651 AbcMatrix::scale(int numberAlreadyScaled)
    652 {
    653   int numberRows=model_->numberRows();
    654   bool inBranchAndBound=(model_->specialOptions(),0x1000000)!=0;
    655   bool doScaling=numberAlreadyScaled>=0;
     643void AbcMatrix::scale(int numberAlreadyScaled)
     644{
     645  int numberRows = model_->numberRows();
     646  bool inBranchAndBound = (model_->specialOptions(), 0x1000000) != 0;
     647  bool doScaling = numberAlreadyScaled >= 0;
    656648  if (!doScaling)
    657     numberAlreadyScaled=0;
    658   if (numberAlreadyScaled==numberRows)
     649    numberAlreadyScaled = 0;
     650  if (numberAlreadyScaled == numberRows)
    659651    return; // no need to do anything
    660   int numberColumns=model_->numberColumns();
    661   double *  COIN_RESTRICT rowScale=model_->rowScale2();
    662   double *  COIN_RESTRICT inverseRowScale=model_->inverseRowScale2();
    663   double *  COIN_RESTRICT columnScale=model_->columnScale2();
    664   double *  COIN_RESTRICT inverseColumnScale=model_->inverseColumnScale2();
     652  int numberColumns = model_->numberColumns();
     653  double *COIN_RESTRICT rowScale = model_->rowScale2();
     654  double *COIN_RESTRICT inverseRowScale = model_->inverseRowScale2();
     655  double *COIN_RESTRICT columnScale = model_->columnScale2();
     656  double *COIN_RESTRICT inverseColumnScale = model_->inverseColumnScale2();
    665657  // we are going to mark bits we are interested in
    666   int whichArray=model_->getAvailableArrayPublic();
    667   char *  COIN_RESTRICT usefulColumn = reinterpret_cast<char *>(model_->usefulArray(whichArray)->getIndices());
    668   memset(usefulColumn,1,numberColumns);
    669   const double *  COIN_RESTRICT rowLower = model_->rowLower();
    670   const double *  COIN_RESTRICT rowUpper = model_->rowUpper();
    671   const double *  COIN_RESTRICT columnLower = model_->columnLower();
    672   const double *  COIN_RESTRICT columnUpper = model_->columnUpper();
     658  int whichArray = model_->getAvailableArrayPublic();
     659  char *COIN_RESTRICT usefulColumn = reinterpret_cast< char * >(model_->usefulArray(whichArray)->getIndices());
     660  memset(usefulColumn, 1, numberColumns);
     661  const double *COIN_RESTRICT rowLower = model_->rowLower();
     662  const double *COIN_RESTRICT rowUpper = model_->rowUpper();
     663  const double *COIN_RESTRICT columnLower = model_->columnLower();
     664  const double *COIN_RESTRICT columnUpper = model_->columnUpper();
    673665  //#define LEAVE_FIXED
    674666  // mark empty and fixed columns
    675667  // get matrix data pointers
    676   const int *  COIN_RESTRICT row = matrix_->getIndices();
    677   const CoinBigIndex *  COIN_RESTRICT columnStart = matrix_->getVectorStarts();
    678   double *  COIN_RESTRICT elementByColumn = matrix_->getMutableElements();
    679   CoinPackedMatrix *  COIN_RESTRICT rowCopy = reverseOrderedCopy();
    680   const int * column = rowCopy->getIndices();
    681   const CoinBigIndex *  COIN_RESTRICT rowStart = rowCopy->getVectorStarts();
    682   const double *  COIN_RESTRICT element = rowCopy->getElements();
    683   assert (numberAlreadyScaled>=0&&numberAlreadyScaled<numberRows);
     668  const int *COIN_RESTRICT row = matrix_->getIndices();
     669  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     670  double *COIN_RESTRICT elementByColumn = matrix_->getMutableElements();
     671  CoinPackedMatrix *COIN_RESTRICT rowCopy = reverseOrderedCopy();
     672  const int *column = rowCopy->getIndices();
     673  const CoinBigIndex *COIN_RESTRICT rowStart = rowCopy->getVectorStarts();
     674  const double *COIN_RESTRICT element = rowCopy->getElements();
     675  assert(numberAlreadyScaled >= 0 && numberAlreadyScaled < numberRows);
    684676#ifndef LEAVE_FIXED
    685677  for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
    686     if (columnUpper[iColumn] == columnLower[iColumn])
    687       usefulColumn[iColumn]=0;
     678    if (columnUpper[iColumn] == columnLower[iColumn])
     679      usefulColumn[iColumn] = 0;
    688680  }
    689681#endif
     
    692684  if (!numberAlreadyScaled) {
    693685    // may be redundant
    694     CoinFillN(rowScale,numberRows,1.0);
    695     CoinFillN(columnScale,numberColumns,1.0);
     686    CoinFillN(rowScale, numberRows, 1.0);
     687    CoinFillN(columnScale, numberColumns, 1.0);
    696688    CoinBigIndex start = 0;
    697689    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
    698       CoinBigIndex end = columnStart[iColumn+1];
     690      CoinBigIndex end = columnStart[iColumn + 1];
    699691      if (usefulColumn[iColumn]) {
    700         if (end>start) {
    701           for (CoinBigIndex j = start; j < end; j++) {
    702             double value = fabs(elementByColumn[j]);
    703             overallLargest = CoinMax(overallLargest, value);
    704             overallSmallest = CoinMin(overallSmallest, value);
    705           }
    706         } else {
    707           usefulColumn[iColumn]=0;
    708         }
    709       }
    710       start=end;
     692        if (end > start) {
     693          for (CoinBigIndex j = start; j < end; j++) {
     694            double value = fabs(elementByColumn[j]);
     695            overallLargest = CoinMax(overallLargest, value);
     696            overallSmallest = CoinMin(overallSmallest, value);
     697          }
     698        } else {
     699          usefulColumn[iColumn] = 0;
     700        }
     701      }
     702      start = end;
    711703    }
    712704  } else {
    713705    CoinBigIndex start = rowStart[numberAlreadyScaled];
    714706    for (int iRow = numberAlreadyScaled; iRow < numberRows; iRow++) {
    715       rowScale[iRow]=1.0;
    716       CoinBigIndex end = rowStart[iRow+1];
     707      rowScale[iRow] = 1.0;
     708      CoinBigIndex end = rowStart[iRow + 1];
    717709      for (CoinBigIndex j = start; j < end; j++) {
    718         int iColumn=column[j];
    719         if (usefulColumn[iColumn]) {
    720           double value = fabs(elementByColumn[j])*columnScale[iColumn];
    721           overallLargest = CoinMax(overallLargest, value);
    722           overallSmallest = CoinMin(overallSmallest, value);
    723         }
    724       }
    725     }
    726   }
    727   if ((overallSmallest >= 0.5 && overallLargest <= 2.0)||!doScaling) {
     710        int iColumn = column[j];
     711        if (usefulColumn[iColumn]) {
     712          double value = fabs(elementByColumn[j]) * columnScale[iColumn];
     713          overallLargest = CoinMax(overallLargest, value);
     714          overallSmallest = CoinMin(overallSmallest, value);
     715        }
     716      }
     717    }
     718  }
     719  if ((overallSmallest >= 0.5 && overallLargest <= 2.0) || !doScaling) {
    728720    //printf("no scaling\n");
    729721    delete rowCopy;
    730722    model_->clearArraysPublic(whichArray);
    731     CoinFillN(inverseRowScale+numberAlreadyScaled,numberRows-numberAlreadyScaled,1.0);
    732     if (!numberAlreadyScaled) 
    733       CoinFillN(inverseColumnScale,numberColumns,1.0);
     723    CoinFillN(inverseRowScale + numberAlreadyScaled, numberRows - numberAlreadyScaled, 1.0);
     724    if (!numberAlreadyScaled)
     725      CoinFillN(inverseColumnScale, numberColumns, 1.0);
    734726    //moveLargestToStart();
    735     return ;
     727    return;
    736728  }
    737729  // need to scale
    738730  double largest;
    739731  double smallest;
    740   int scalingMethod=model_->scalingFlag();
     732  int scalingMethod = model_->scalingFlag();
    741733  if (scalingMethod == 4) {
    742734    // As auto
     
    751743  // if scalingMethod 3 then may change
    752744  bool extraDetails = (model_->logLevel() > 2);
    753   bool secondTime=false;
     745  bool secondTime = false;
    754746  while (!finished) {
    755747    int numberPass = !numberAlreadyScaled ? 3 : 1;
     
    757749    overallSmallest = 1.0e20;
    758750    if (!secondTime) {
    759       secondTime=true;
     751      secondTime = true;
    760752    } else {
    761       CoinFillN ( rowScale, numberRows, 1.0);
    762       CoinFillN ( columnScale, numberColumns, 1.0);
    763     } 
     753      CoinFillN(rowScale, numberRows, 1.0);
     754      CoinFillN(columnScale, numberColumns, 1.0);
     755    }
    764756    if (scalingMethod == 1 || scalingMethod == 3) {
    765757      // Maximum in each row
    766758      for (int iRow = numberAlreadyScaled; iRow < numberRows; iRow++) {
    767         largest = 1.0e-10;
    768         for (CoinBigIndex j = rowStart[iRow]; j < rowStart[iRow+1]; j++) {
    769           int iColumn = column[j];
    770           if (usefulColumn[iColumn]) {
    771             double value = fabs(element[j]);
    772             largest = CoinMax(largest, value);
    773             assert (largest < 1.0e40);
    774           }
    775         }
    776         rowScale[iRow] = 1.0 / largest;
     759        largest = 1.0e-10;
     760        for (CoinBigIndex j = rowStart[iRow]; j < rowStart[iRow + 1]; j++) {
     761          int iColumn = column[j];
     762          if (usefulColumn[iColumn]) {
     763            double value = fabs(element[j]);
     764            largest = CoinMax(largest, value);
     765            assert(largest < 1.0e40);
     766          }
     767        }
     768        rowScale[iRow] = 1.0 / largest;
    777769#ifdef COIN_DEVELOP
    778         if (extraDetails) {
    779           overallLargest = CoinMax(overallLargest, largest);
    780           overallSmallest = CoinMin(overallSmallest, largest);
    781         }
     770        if (extraDetails) {
     771          overallLargest = CoinMax(overallLargest, largest);
     772          overallSmallest = CoinMin(overallSmallest, largest);
     773        }
    782774#endif
    783775      }
    784776    } else {
    785777      while (numberPass) {
    786         overallLargest = 0.0;
    787         overallSmallest = 1.0e50;
    788         numberPass--;
    789         // Geometric mean on row scales
    790         for (int iRow = numberAlreadyScaled; iRow < numberRows; iRow++) {
    791           largest = 1.0e-50;
    792           smallest = 1.0e50;
    793           for (CoinBigIndex j = rowStart[iRow]; j < rowStart[iRow+1]; j++) {
    794             int iColumn = column[j];
    795             if (usefulColumn[iColumn]) {
    796               double value = fabs(element[j]);
    797               value *= columnScale[iColumn];
    798               largest = CoinMax(largest, value);
    799               smallest = CoinMin(smallest, value);
    800             }
    801           }
     778        overallLargest = 0.0;
     779        overallSmallest = 1.0e50;
     780        numberPass--;
     781        // Geometric mean on row scales
     782        for (int iRow = numberAlreadyScaled; iRow < numberRows; iRow++) {
     783          largest = 1.0e-50;
     784          smallest = 1.0e50;
     785          for (CoinBigIndex j = rowStart[iRow]; j < rowStart[iRow + 1]; j++) {
     786            int iColumn = column[j];
     787            if (usefulColumn[iColumn]) {
     788              double value = fabs(element[j]);
     789              value *= columnScale[iColumn];
     790              largest = CoinMax(largest, value);
     791              smallest = CoinMin(smallest, value);
     792            }
     793          }
    802794#ifdef SQRT_ARRAY
    803           rowScale[iRow] = smallest * largest;
     795          rowScale[iRow] = smallest * largest;
    804796#else
    805           rowScale[iRow] = 1.0 / sqrt(smallest * largest);
    806 #endif
    807           if (extraDetails) {
    808             overallLargest = CoinMax(largest * rowScale[iRow], overallLargest);
    809             overallSmallest = CoinMin(smallest * rowScale[iRow], overallSmallest);
    810           }
    811         }
    812         if (model_->scalingFlag() == 5)
    813           break; // just scale rows
     797          rowScale[iRow] = 1.0 / sqrt(smallest * largest);
     798#endif
     799          if (extraDetails) {
     800            overallLargest = CoinMax(largest * rowScale[iRow], overallLargest);
     801            overallSmallest = CoinMin(smallest * rowScale[iRow], overallSmallest);
     802          }
     803        }
     804        if (model_->scalingFlag() == 5)
     805          break; // just scale rows
    814806#ifdef SQRT_ARRAY
    815         CoinAbcInverseSqrts(rowScale, numberRows);
    816 #endif
    817         if (!inBranchAndBound)
    818         model_->messageHandler()->message(CLP_PACKEDSCALE_WHILE, *model_->messagesPointer())
    819           << overallSmallest
    820           << overallLargest
    821           << CoinMessageEol;
    822         // skip last column round
    823         if (numberPass == 1)
    824           break;
    825         // Geometric mean on column scales
    826         for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
    827           if (usefulColumn[iColumn]) {
    828             largest = 1.0e-50;
    829             smallest = 1.0e50;
    830             for (CoinBigIndex j = columnStart[iColumn];
    831                  j < columnStart[iColumn+1]; j++) {
    832               int iRow = row[j];
    833               double value = fabs(elementByColumn[j]);
    834               value *= rowScale[iRow];
    835               largest = CoinMax(largest, value);
    836               smallest = CoinMin(smallest, value);
    837             }
     807        CoinAbcInverseSqrts(rowScale, numberRows);
     808#endif
     809        if (!inBranchAndBound)
     810          model_->messageHandler()->message(CLP_PACKEDSCALE_WHILE, *model_->messagesPointer())
     811            << overallSmallest
     812            << overallLargest
     813            << CoinMessageEol;
     814        // skip last column round
     815        if (numberPass == 1)
     816          break;
     817        // Geometric mean on column scales
     818        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     819          if (usefulColumn[iColumn]) {
     820            largest = 1.0e-50;
     821            smallest = 1.0e50;
     822            for (CoinBigIndex j = columnStart[iColumn];
     823                 j < columnStart[iColumn + 1]; j++) {
     824              int iRow = row[j];
     825              double value = fabs(elementByColumn[j]);
     826              value *= rowScale[iRow];
     827              largest = CoinMax(largest, value);
     828              smallest = CoinMin(smallest, value);
     829            }
    838830#ifdef SQRT_ARRAY
    839             columnScale[iColumn] = smallest * largest;
     831            columnScale[iColumn] = smallest * largest;
    840832#else
    841             columnScale[iColumn] = 1.0 / sqrt(smallest * largest);
    842 #endif
    843           }
    844         }
     833            columnScale[iColumn] = 1.0 / sqrt(smallest * largest);
     834#endif
     835          }
     836        }
    845837#ifdef SQRT_ARRAY
    846         CoinAbcInverseSqrts(columnScale, numberColumns);
     838        CoinAbcInverseSqrts(columnScale, numberColumns);
    847839#endif
    848840      }
     
    853845      double scaledDifference = difference * rowScale[iRow];
    854846      if (scaledDifference > tolerance && scaledDifference < 1.0e-4) {
    855         // make gap larger
    856         rowScale[iRow] *= 1.0e-4 / scaledDifference;
    857         rowScale[iRow] = CoinMax(1.0e-10, CoinMin(1.0e10, rowScale[iRow]));
    858         //printf("Row %d difference %g scaled diff %g => %g\n",iRow,difference,
    859         // scaledDifference,difference*rowScale[iRow]);
     847        // make gap larger
     848        rowScale[iRow] *= 1.0e-4 / scaledDifference;
     849        rowScale[iRow] = CoinMax(1.0e-10, CoinMin(1.0e10, rowScale[iRow]));
     850        //printf("Row %d difference %g scaled diff %g => %g\n",iRow,difference,
     851        // scaledDifference,difference*rowScale[iRow]);
    860852      }
    861853    }
     
    865857      overallSmallest = 1.0e50;
    866858      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
    867         if (usefulColumn[iColumn]) {
    868           largest = 1.0e-20;
    869           smallest = 1.0e50;
    870           for (CoinBigIndex j = columnStart[iColumn];j < columnStart[iColumn+1]; j++) {
    871             int iRow = row[j];
    872             double value = fabs(elementByColumn[j] * rowScale[iRow]);
    873             largest = CoinMax(largest, value);
    874             smallest = CoinMin(smallest, value);
    875           }
    876           if (overallSmallest * largest > smallest)
    877             overallSmallest = smallest / largest;
    878         }
     859        if (usefulColumn[iColumn]) {
     860          largest = 1.0e-20;
     861          smallest = 1.0e50;
     862          for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn + 1]; j++) {
     863            int iRow = row[j];
     864            double value = fabs(elementByColumn[j] * rowScale[iRow]);
     865            largest = CoinMax(largest, value);
     866            smallest = CoinMin(smallest, value);
     867          }
     868          if (overallSmallest * largest > smallest)
     869            overallSmallest = smallest / largest;
     870        }
    879871      }
    880872    }
     
    885877      scalingMethod = 4;
    886878    } else {
    887       assert (scalingMethod == 4);
     879      assert(scalingMethod == 4);
    888880      if (overallSmallest > 2.0 * savedOverallRatio) {
    889         finished = true; // geometric was better
    890         if (model_->scalingFlag() == 4) {
    891           // If in Branch and bound change
    892           if ((model_->specialOptions() & 1024) != 0) {
    893             model_->scaling(2);
    894           }
    895         }
     881        finished = true; // geometric was better
     882        if (model_->scalingFlag() == 4) {
     883          // If in Branch and bound change
     884          if ((model_->specialOptions() & 1024) != 0) {
     885            model_->scaling(2);
     886          }
     887        }
    896888      } else {
    897         scalingMethod = 1; // redo equilibrium
    898         if (model_->scalingFlag() == 4) {
    899           // If in Branch and bound change
    900           if ((model_->specialOptions() & 1024) != 0) {
    901             model_->scaling(1);
    902           }
    903         }
     889        scalingMethod = 1; // redo equilibrium
     890        if (model_->scalingFlag() == 4) {
     891          // If in Branch and bound change
     892          if ((model_->specialOptions() & 1024) != 0) {
     893            model_->scaling(1);
     894          }
     895        }
    904896      }
    905897#if 0
     
    928920  overallLargest = CoinMin(100.0, overallLargest);
    929921  overallSmallest = 1.0e50;
    930   char *  COIN_RESTRICT usedRow = reinterpret_cast<char *>(inverseRowScale);
     922  char *COIN_RESTRICT usedRow = reinterpret_cast< char * >(inverseRowScale);
    931923  memset(usedRow, 0, numberRows);
    932924  //printf("scaling %d\n",model_->scalingFlag());
     
    934926    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
    935927      if (usefulColumn[iColumn]) {
    936         largest = 1.0e-20;
    937         smallest = 1.0e50;
    938         for (CoinBigIndex j = columnStart[iColumn];j < columnStart[iColumn+1]; j++) {
    939           int iRow = row[j];
    940           usedRow[iRow] = 1;
    941           double value = fabs(elementByColumn[j] * rowScale[iRow]);
    942           largest = CoinMax(largest, value);
    943           smallest = CoinMin(smallest, value);
    944         }
    945         columnScale[iColumn] = overallLargest / largest;
    946         //columnScale[iColumn]=CoinMax(1.0e-10,CoinMin(1.0e10,columnScale[iColumn]));
     928        largest = 1.0e-20;
     929        smallest = 1.0e50;
     930        for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn + 1]; j++) {
     931          int iRow = row[j];
     932          usedRow[iRow] = 1;
     933          double value = fabs(elementByColumn[j] * rowScale[iRow]);
     934          largest = CoinMax(largest, value);
     935          smallest = CoinMin(smallest, value);
     936        }
     937        columnScale[iColumn] = overallLargest / largest;
     938        //columnScale[iColumn]=CoinMax(1.0e-10,CoinMin(1.0e10,columnScale[iColumn]));
    947939#ifdef RANDOMIZE
    948         if (!numberAlreadyScaled) {
    949           double value = 0.5 - randomNumberGenerator_.randomDouble(); //between -0.5 to + 0.5
    950           columnScale[iColumn] *= (1.0 + 0.1 * value);
    951         }
    952 #endif
    953         double difference = columnUpper[iColumn] - columnLower[iColumn];
    954         if (difference < 1.0e-5 * columnScale[iColumn]) {
    955           // make gap larger
    956           columnScale[iColumn] = difference / 1.0e-5;
    957           //printf("Column %d difference %g scaled diff %g => %g\n",iColumn,difference,
    958           // scaledDifference,difference*columnScale[iColumn]);
    959         }
    960         double value = smallest * columnScale[iColumn];
    961         if (overallSmallest > value)
    962           overallSmallest = value;
    963         //overallSmallest = CoinMin(overallSmallest,smallest*columnScale[iColumn]);
     940        if (!numberAlreadyScaled) {
     941          double value = 0.5 - randomNumberGenerator_.randomDouble(); //between -0.5 to + 0.5
     942          columnScale[iColumn] *= (1.0 + 0.1 * value);
     943        }
     944#endif
     945        double difference = columnUpper[iColumn] - columnLower[iColumn];
     946        if (difference < 1.0e-5 * columnScale[iColumn]) {
     947          // make gap larger
     948          columnScale[iColumn] = difference / 1.0e-5;
     949          //printf("Column %d difference %g scaled diff %g => %g\n",iColumn,difference,
     950          // scaledDifference,difference*columnScale[iColumn]);
     951        }
     952        double value = smallest * columnScale[iColumn];
     953        if (overallSmallest > value)
     954          overallSmallest = value;
     955        //overallSmallest = CoinMin(overallSmallest,smallest*columnScale[iColumn]);
    964956      } else {
    965         assert(columnScale[iColumn] == 1.0);
    966         //columnScale[iColumn]=1.0;
     957        assert(columnScale[iColumn] == 1.0);
     958        //columnScale[iColumn]=1.0;
    967959      }
    968960    }
    969961    for (int iRow = numberAlreadyScaled; iRow < numberRows; iRow++) {
    970962      if (!usedRow[iRow]) {
    971         rowScale[iRow] = 1.0;
     963        rowScale[iRow] = 1.0;
    972964      }
    973965    }
    974966  }
    975967  if (!inBranchAndBound)
    976   model_->messageHandler()->message(CLP_PACKEDSCALE_FINAL, *model_->messagesPointer())
    977     << overallSmallest
    978     << overallLargest
    979     << CoinMessageEol;
     968    model_->messageHandler()->message(CLP_PACKEDSCALE_FINAL, *model_->messagesPointer())
     969      << overallSmallest
     970      << overallLargest
     971      << CoinMessageEol;
    980972  if (overallSmallest < 1.0e-13) {
    981973    // Change factorization zero tolerance
    982974    double newTolerance = CoinMax(1.0e-15 * (overallSmallest / 1.0e-13),
    983                                   1.0e-18);
     975      1.0e-18);
    984976    if (model_->factorization()->zeroTolerance() > newTolerance)
    985977      model_->factorization()->zeroTolerance(newTolerance);
     
    987979    model_->setZeroTolerance(newTolerance);
    988980#ifndef NDEBUG
    989     assert (newTolerance<0.0); // just so we can fix
     981    assert(newTolerance < 0.0); // just so we can fix
    990982#endif
    991983  }
    992984  // make copy (could do faster by using previous values)
    993985  // could just do partial
    994   CoinAbcReciprocal(inverseRowScale+numberAlreadyScaled,numberRows-numberAlreadyScaled,
    995                 rowScale+numberAlreadyScaled);
    996   if (!numberAlreadyScaled) 
    997     CoinAbcReciprocal(inverseColumnScale,numberColumns,columnScale);
     986  CoinAbcReciprocal(inverseRowScale + numberAlreadyScaled, numberRows - numberAlreadyScaled,
     987    rowScale + numberAlreadyScaled);
     988  if (!numberAlreadyScaled)
     989    CoinAbcReciprocal(inverseColumnScale, numberColumns, columnScale);
    998990  // Do scaled copy //NO and move largest to start
    999991  CoinBigIndex start = 0;
    1000992  for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
    1001     double scale=columnScale[iColumn];
    1002     CoinBigIndex end = columnStart[iColumn+1];
     993    double scale = columnScale[iColumn];
     994    CoinBigIndex end = columnStart[iColumn + 1];
    1003995    for (CoinBigIndex j = start; j < end; j++) {
    1004       double value=elementByColumn[j];
     996      double value = elementByColumn[j];
    1005997      int iRow = row[j];
    1006       if (iRow>=numberAlreadyScaled) {
    1007         value*= scale*rowScale[iRow];
    1008         elementByColumn[j]= value;
    1009       }
    1010     }
    1011     start=end;
     998      if (iRow >= numberAlreadyScaled) {
     999        value *= scale * rowScale[iRow];
     1000        elementByColumn[j] = value;
     1001      }
     1002    }
     1003    start = end;
    10121004  }
    10131005  delete rowCopy;
     
    10411033   @pre <code>y</code> must be of size <code>numRows()</code> */
    10421034//scaled versions
    1043 void
    1044 AbcMatrix::timesModifyExcludingSlacks(double scalar,
    1045                                       const double * x, double * y) const
     1035void AbcMatrix::timesModifyExcludingSlacks(double scalar,
     1036  const double *x, double *y) const
    10461037{
    10471038  int numberTotal = model_->numberTotal();
    10481039  int maximumRows = model_->maximumAbcNumberRows();
    10491040  // get matrix data pointers
    1050   const int *  COIN_RESTRICT row = matrix_->getIndices();
    1051   const CoinBigIndex *  COIN_RESTRICT columnStart = matrix_->getVectorStarts()-maximumRows;
    1052   const double *  COIN_RESTRICT elementByColumn = matrix_->getElements();
     1041  const int *COIN_RESTRICT row = matrix_->getIndices();
     1042  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts() - maximumRows;
     1043  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
    10531044  for (int iColumn = maximumRows; iColumn < numberTotal; iColumn++) {
    10541045    double value = x[iColumn];
    10551046    if (value) {
    10561047      CoinBigIndex start = columnStart[iColumn];
    1057       CoinBigIndex end = columnStart[iColumn+1];
     1048      CoinBigIndex end = columnStart[iColumn + 1];
    10581049      value *= scalar;
    10591050      for (CoinBigIndex j = start; j < end; j++) {
    1060         int iRow = row[j];
    1061         y[iRow] += value * elementByColumn[j];
     1051        int iRow = row[j];
     1052        y[iRow] += value * elementByColumn[j];
    10621053      }
    10631054    }
     
    10671058   @pre <code>x</code> must be of size <code>numColumns()+numRows()</code>
    10681059   @pre <code>y</code> must be of size <code>numRows()</code> */
    1069 void
    1070 AbcMatrix::timesModifyIncludingSlacks(double scalar,
    1071                                       const double * x, double * y) const
    1072 {
    1073   int numberRows=model_->numberRows();
     1060void AbcMatrix::timesModifyIncludingSlacks(double scalar,
     1061  const double *x, double *y) const
     1062{
     1063  int numberRows = model_->numberRows();
    10741064  int numberTotal = model_->numberTotal();
    10751065  int maximumRows = model_->maximumAbcNumberRows();
    10761066  // makes no sense for x==y??
    1077   assert (x!=y);
     1067  assert(x != y);
    10781068  // For now just by column and assumes already scaled (use reallyScale)
    10791069  // get matrix data pointers
    1080   const int * COIN_RESTRICT row = matrix_->getIndices();
    1081   const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts()-maximumRows;
    1082   const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
    1083   if (scalar==1.0) {
     1070  const int *COIN_RESTRICT row = matrix_->getIndices();
     1071  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts() - maximumRows;
     1072  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
     1073  if (scalar == 1.0) {
    10841074    // add
    1085     if (x!=y) {
    1086       for (int i=0;i<numberRows;i++)
    1087         y[i]+=x[i];
     1075    if (x != y) {
     1076      for (int i = 0; i < numberRows; i++)
     1077        y[i] += x[i];
    10881078    } else {
    1089       for (int i=0;i<numberRows;i++)
    1090         y[i]+=y[i];
     1079      for (int i = 0; i < numberRows; i++)
     1080        y[i] += y[i];
    10911081    }
    10921082    for (int iColumn = maximumRows; iColumn < numberTotal; iColumn++) {
    10931083      double value = x[iColumn];
    10941084      if (value) {
    1095         CoinBigIndex start = columnStart[iColumn];
    1096         CoinBigIndex next = columnStart[iColumn+1];
    1097         for (CoinBigIndex j = start; j < next; j++) {
    1098           int jRow = row[j];
    1099           y[jRow]+=value*elementByColumn[j];
    1100         }
     1085        CoinBigIndex start = columnStart[iColumn];
     1086        CoinBigIndex next = columnStart[iColumn + 1];
     1087        for (CoinBigIndex j = start; j < next; j++) {
     1088          int jRow = row[j];
     1089          y[jRow] += value * elementByColumn[j];
     1090        }
    11011091      }
    11021092    }
    11031093  } else {
    11041094    // subtract
    1105     if (x!=y) {
    1106       for (int i=0;i<numberRows;i++)
    1107         y[i]-=x[i];
     1095    if (x != y) {
     1096      for (int i = 0; i < numberRows; i++)
     1097        y[i] -= x[i];
    11081098    } else {
    1109       for (int i=0;i<numberRows;i++)
    1110         y[i]=0.0;
     1099      for (int i = 0; i < numberRows; i++)
     1100        y[i] = 0.0;
    11111101    }
    11121102    for (int iColumn = maximumRows; iColumn < numberTotal; iColumn++) {
    11131103      double value = x[iColumn];
    11141104      if (value) {
    1115         CoinBigIndex start = columnStart[iColumn];
    1116         CoinBigIndex next = columnStart[iColumn+1];
    1117         for (CoinBigIndex j = start; j < next; j++) {
    1118           int jRow = row[j];
    1119           y[jRow]-=value*elementByColumn[j];
    1120         }
     1105        CoinBigIndex start = columnStart[iColumn];
     1106        CoinBigIndex next = columnStart[iColumn + 1];
     1107        for (CoinBigIndex j = start; j < next; j++) {
     1108          int jRow = row[j];
     1109          y[jRow] -= value * elementByColumn[j];
     1110        }
    11211111      }
    11221112    }
     
    11261116   @pre <code>x</code> must be of size <code>numColumns()+numRows()</code>
    11271117   @pre <code>y</code> must be of size <code>numRows()</code> */
    1128 void
    1129 AbcMatrix::timesIncludingSlacks(double scalar,
    1130                                 const double * x, double * y) const
    1131 {
    1132   int numberRows=model_->numberRows();
     1118void AbcMatrix::timesIncludingSlacks(double scalar,
     1119  const double *x, double *y) const
     1120{
     1121  int numberRows = model_->numberRows();
    11331122  int numberTotal = model_->numberTotal();
    11341123  int maximumRows = model_->maximumAbcNumberRows();
    11351124  // For now just by column and assumes already scaled (use reallyScale)
    11361125  // get matrix data pointers
    1137   const int * COIN_RESTRICT row = matrix_->getIndices();
    1138   const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts()-maximumRows;
    1139   const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
    1140   if (scalar==1.0) {
     1126  const int *COIN_RESTRICT row = matrix_->getIndices();
     1127  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts() - maximumRows;
     1128  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
     1129  if (scalar == 1.0) {
    11411130    // add
    1142     if (x!=y) {
    1143       for (int i=0;i<numberRows;i++)
    1144         y[i]=x[i];
    1145     } 
     1131    if (x != y) {
     1132      for (int i = 0; i < numberRows; i++)
     1133        y[i] = x[i];
     1134    }
    11461135    for (int iColumn = maximumRows; iColumn < numberTotal; iColumn++) {
    11471136      double value = x[iColumn];
    11481137      if (value) {
    1149         CoinBigIndex start = columnStart[iColumn];
    1150         CoinBigIndex next = columnStart[iColumn+1];
    1151         for (CoinBigIndex j = start; j < next; j++) {
    1152           int jRow = row[j];
    1153           y[jRow]+=value*elementByColumn[j];
    1154         }
     1138        CoinBigIndex start = columnStart[iColumn];
     1139        CoinBigIndex next = columnStart[iColumn + 1];
     1140        for (CoinBigIndex j = start; j < next; j++) {
     1141          int jRow = row[j];
     1142          y[jRow] += value * elementByColumn[j];
     1143        }
    11551144      }
    11561145    }
    11571146  } else {
    11581147    // subtract
    1159     if (x!=y) {
    1160       for (int i=0;i<numberRows;i++)
    1161         y[i]=-x[i];
     1148    if (x != y) {
     1149      for (int i = 0; i < numberRows; i++)
     1150        y[i] = -x[i];
    11621151    } else {
    1163       for (int i=0;i<numberRows;i++)
    1164         y[i]=-y[i];
     1152      for (int i = 0; i < numberRows; i++)
     1153        y[i] = -y[i];
    11651154    }
    11661155    for (int iColumn = maximumRows; iColumn < numberTotal; iColumn++) {
    11671156      double value = x[iColumn];
    11681157      if (value) {
    1169         CoinBigIndex start = columnStart[iColumn];
    1170         CoinBigIndex next = columnStart[iColumn+1];
    1171         for (CoinBigIndex j = start; j < next; j++) {
    1172           int jRow = row[j];
    1173           y[jRow]-=value*elementByColumn[j];
    1174         }
    1175       }
    1176     }
    1177   }
    1178 }
    1179 extern double parallelDual4(AbcSimplexDual * );
    1180 static double firstPass(AbcSimplex * model,int iBlock,
    1181                  double upperTheta, int & freeSequence,
    1182                  CoinPartitionedVector & tableauRow,
    1183                  CoinPartitionedVector & candidateList)
    1184 {
    1185   int *  COIN_RESTRICT index = tableauRow.getIndices();
    1186   double *  COIN_RESTRICT array = tableauRow.denseVector();
    1187   double *  COIN_RESTRICT arrayCandidate=candidateList.denseVector();
    1188   int *  COIN_RESTRICT indexCandidate = candidateList.getIndices();
    1189   const double *  COIN_RESTRICT abcDj = model->djRegion();
    1190   const unsigned char *  COIN_RESTRICT internalStatus = model->internalStatus();
     1158        CoinBigIndex start = columnStart[iColumn];
     1159        CoinBigIndex next = columnStart[iColumn + 1];
     1160        for (CoinBigIndex j = start; j < next; j++) {
     1161          int jRow = row[j];
     1162          y[jRow] -= value * elementByColumn[j];
     1163        }
     1164      }
     1165    }
     1166  }
     1167}
     1168extern double parallelDual4(AbcSimplexDual *);
     1169static double firstPass(AbcSimplex *model, int iBlock,
     1170  double upperTheta, int &freeSequence,
     1171  CoinPartitionedVector &tableauRow,
     1172  CoinPartitionedVector &candidateList)
     1173{
     1174  int *COIN_RESTRICT index = tableauRow.getIndices();
     1175  double *COIN_RESTRICT array = tableauRow.denseVector();
     1176  double *COIN_RESTRICT arrayCandidate = candidateList.denseVector();
     1177  int *COIN_RESTRICT indexCandidate = candidateList.getIndices();
     1178  const double *COIN_RESTRICT abcDj = model->djRegion();
     1179  const unsigned char *COIN_RESTRICT internalStatus = model->internalStatus();
    11911180  // do first pass to get possibles
    11921181  double bestPossible = 0.0;
     
    11941183  double tentativeTheta = 1.0e25; // try with this much smaller as guess
    11951184  double acceptablePivot = model->currentAcceptablePivot();
    1196   double dualT=-model->currentDualTolerance();
     1185  double dualT = -model->currentDualTolerance();
    11971186  // fixed will have been taken out by now
    1198   const double multiplier[] = { 1.0, -1.0};
    1199   freeSequence=-1;
    1200   int firstIn=model->abcMatrix()->blockStart(iBlock);
    1201   int numberNonZero=tableauRow.getNumElements(iBlock)+firstIn;
    1202   int numberRemaining=firstIn;
     1187  const double multiplier[] = { 1.0, -1.0 };
     1188  freeSequence = -1;
     1189  int firstIn = model->abcMatrix()->blockStart(iBlock);
     1190  int numberNonZero = tableauRow.getNumElements(iBlock) + firstIn;
     1191  int numberRemaining = firstIn;
    12031192  //first=tableauRow.getNumElements();
    12041193  // could pass in last and if numberNonZero==last-firstIn scan as well
     
    12061195    for (int i = firstIn; i < numberNonZero; i++) {
    12071196      int iSequence = index[i];
    1208       double tableauValue=array[i];
    1209       unsigned char iStatus=internalStatus[iSequence]&7;
     1197      double tableauValue = array[i];
     1198      unsigned char iStatus = internalStatus[iSequence] & 7;
    12101199      double mult = multiplier[iStatus];
    12111200      double alpha = tableauValue * mult;
     
    12131202      double value = oldValue - tentativeTheta * alpha;
    12141203      if (value < dualT) {
    1215         bestPossible = CoinMax(bestPossible, alpha);
    1216         value = oldValue - upperTheta * alpha;
    1217         if (value < dualT && alpha >= acceptablePivot) {
    1218           upperTheta = (oldValue - dualT) / alpha;
    1219         }
    1220         // add to list
    1221         arrayCandidate[numberRemaining] = alpha;
    1222         indexCandidate[numberRemaining++] = iSequence;
     1204        bestPossible = CoinMax(bestPossible, alpha);
     1205        value = oldValue - upperTheta * alpha;
     1206        if (value < dualT && alpha >= acceptablePivot) {
     1207          upperTheta = (oldValue - dualT) / alpha;
     1208        }
     1209        // add to list
     1210        arrayCandidate[numberRemaining] = alpha;
     1211        indexCandidate[numberRemaining++] = iSequence;
    12231212      }
    12241213    }
     
    12261215    double badFree = 0.0;
    12271216    double freeAlpha = model->currentAcceptablePivot();
    1228     int freeSequenceIn=model->freeSequenceIn();
     1217    int freeSequenceIn = model->freeSequenceIn();
    12291218    double currentDualTolerance = model->currentDualTolerance();
    12301219    for (int i = firstIn; i < numberNonZero; i++) {
    12311220      int iSequence = index[i];
    1232       double tableauValue=array[i];
    1233       unsigned char iStatus=internalStatus[iSequence]&7;
    1234       if ((iStatus&4)==0) {
    1235         double mult = multiplier[iStatus];
    1236         double alpha = tableauValue * mult;
    1237         double oldValue = abcDj[iSequence] * mult;
    1238         double value = oldValue - tentativeTheta * alpha;
    1239         if (value < dualT) {
    1240           bestPossible = CoinMax(bestPossible, alpha);
    1241           value = oldValue - upperTheta * alpha;
    1242           if (value < dualT && alpha >= acceptablePivot) {
    1243             upperTheta = (oldValue - dualT) / alpha;
    1244           }
    1245           // add to list
    1246           arrayCandidate[numberRemaining] = alpha;
    1247           indexCandidate[numberRemaining++] = iSequence;
    1248         }
     1221      double tableauValue = array[i];
     1222      unsigned char iStatus = internalStatus[iSequence] & 7;
     1223      if ((iStatus & 4) == 0) {
     1224        double mult = multiplier[iStatus];
     1225        double alpha = tableauValue * mult;
     1226        double oldValue = abcDj[iSequence] * mult;
     1227        double value = oldValue - tentativeTheta * alpha;
     1228        if (value < dualT) {
     1229          bestPossible = CoinMax(bestPossible, alpha);
     1230          value = oldValue - upperTheta * alpha;
     1231          if (value < dualT && alpha >= acceptablePivot) {
     1232            upperTheta = (oldValue - dualT) / alpha;
     1233          }
     1234          // add to list
     1235          arrayCandidate[numberRemaining] = alpha;
     1236          indexCandidate[numberRemaining++] = iSequence;
     1237        }
    12491238      } else {
    1250         bool keep;
    1251         bestPossible = CoinMax(bestPossible, fabs(tableauValue));
    1252         double oldValue = abcDj[iSequence];
    1253         // If free has to be very large - should come in via dualRow
    1254         //if (getInternalStatus(iSequence+addSequence)==isFree&&fabs(tableauValue)<1.0e-3)
    1255         //break;
    1256         if (oldValue > currentDualTolerance) {
    1257           keep = true;
    1258         } else if (oldValue < -currentDualTolerance) {
    1259           keep = true;
    1260         } else {
    1261           if (fabs(tableauValue) > CoinMax(10.0 * acceptablePivot, 1.0e-5)) {
    1262             keep = true;
    1263           } else {
    1264             keep = false;
    1265             badFree = CoinMax(badFree, fabs(tableauValue));
    1266           }
    1267         }
    1268         if (keep) {
     1239        bool keep;
     1240        bestPossible = CoinMax(bestPossible, fabs(tableauValue));
     1241        double oldValue = abcDj[iSequence];
     1242        // If free has to be very large - should come in via dualRow
     1243        //if (getInternalStatus(iSequence+addSequence)==isFree&&fabs(tableauValue)<1.0e-3)
     1244        //break;
     1245        if (oldValue > currentDualTolerance) {
     1246          keep = true;
     1247        } else if (oldValue < -currentDualTolerance) {
     1248          keep = true;
     1249        } else {
     1250          if (fabs(tableauValue) > CoinMax(10.0 * acceptablePivot, 1.0e-5)) {
     1251            keep = true;
     1252          } else {
     1253            keep = false;
     1254            badFree = CoinMax(badFree, fabs(tableauValue));
     1255          }
     1256        }
     1257        if (keep) {
    12691258#ifdef PAN
    1270           if (model->fakeSuperBasic(iSequence)>=0) {
    1271 #endif
    1272             if (iSequence==freeSequenceIn)
    1273               tableauValue=COIN_DBL_MAX;
    1274             // free - choose largest
    1275             if (fabs(tableauValue) > fabs(freeAlpha)) {
    1276               freeAlpha = tableauValue;
    1277               freeSequence = iSequence;
    1278             }
     1259          if (model->fakeSuperBasic(iSequence) >= 0) {
     1260#endif
     1261            if (iSequence == freeSequenceIn)
     1262              tableauValue = COIN_DBL_MAX;
     1263            // free - choose largest
     1264            if (fabs(tableauValue) > fabs(freeAlpha)) {
     1265              freeAlpha = tableauValue;
     1266              freeSequence = iSequence;
     1267            }
    12791268#ifdef PAN
    1280           }
    1281 #endif
    1282         }
     1269          }
     1270#endif
     1271        }
    12831272      }
    12841273    }
     
    12861275  //firstInX=numberNonZero-firstIn;
    12871276  //lastInX=-1;//numberRemaining-lastInX;
    1288   tableauRow.setNumElementsPartition(iBlock,numberNonZero-firstIn);
    1289   candidateList.setNumElementsPartition(iBlock,numberRemaining-firstIn);
     1277  tableauRow.setNumElementsPartition(iBlock, numberNonZero - firstIn);
     1278  candidateList.setNumElementsPartition(iBlock, numberRemaining - firstIn);
    12901279  return upperTheta;
    12911280}
    12921281// gets sorted tableau row and a possible value of theta
    1293 double 
    1294 AbcMatrix::dualColumn1Row(int iBlock, double upperTheta, int & freeSequence,
    1295                           const CoinIndexedVector & update,
    1296                           CoinPartitionedVector & tableauRow,
    1297                           CoinPartitionedVector & candidateList) const
     1282double
     1283AbcMatrix::dualColumn1Row(int iBlock, double upperTheta, int &freeSequence,
     1284  const CoinIndexedVector &update,
     1285  CoinPartitionedVector &tableauRow,
     1286  CoinPartitionedVector &candidateList) const
    12981287{
    12991288  int maximumRows = model_->maximumAbcNumberRows();
    1300   int number=update.getNumElements();
    1301   const double *  COIN_RESTRICT pi=update.denseVector();
    1302   const int *  COIN_RESTRICT piIndex = update.getIndices();
    1303   const CoinBigIndex * COIN_RESTRICT rowStart = rowStart_;
    1304   int numberRows=model_->numberRows();
    1305   const CoinBigIndex * COIN_RESTRICT rowEnd = rowStart+numberRows*numberRowBlocks_;
     1289  int number = update.getNumElements();
     1290  const double *COIN_RESTRICT pi = update.denseVector();
     1291  const int *COIN_RESTRICT piIndex = update.getIndices();
     1292  const CoinBigIndex *COIN_RESTRICT rowStart = rowStart_;
     1293  int numberRows = model_->numberRows();
     1294  const CoinBigIndex *COIN_RESTRICT rowEnd = rowStart + numberRows * numberRowBlocks_;
    13061295  // count down
    13071296  int nColumns;
    1308   int firstIn=blockStart_[iBlock];
    1309   int first=firstIn;
     1297  int firstIn = blockStart_[iBlock];
     1298  int first = firstIn;
    13101299  if (!first)
    1311     first=maximumRows;
    1312   int last=blockStart_[iBlock+1];
    1313   nColumns=last-first;
    1314   int target=nColumns;
    1315   rowStart += iBlock*numberRows;
    1316   rowEnd = rowStart+numberRows;
    1317   for (int i=0;i<number;i++) {
    1318     int iRow=piIndex[i];
     1300    first = maximumRows;
     1301  int last = blockStart_[iBlock + 1];
     1302  nColumns = last - first;
     1303  int target = nColumns;
     1304  rowStart += iBlock * numberRows;
     1305  rowEnd = rowStart + numberRows;
     1306  for (int i = 0; i < number; i++) {
     1307    int iRow = piIndex[i];
    13191308    CoinBigIndex end = rowEnd[iRow];
    1320     target -= end-rowStart[iRow];
    1321     if (target<0)
     1309    target -= end - rowStart[iRow];
     1310    if (target < 0)
    13221311      break;
    13231312  }
    1324   if (target>0) {
     1313  if (target > 0) {
    13251314    //printf("going to few %d ops %d\n",number,nColumns-target);
    13261315    return dualColumn1RowFew(iBlock, upperTheta, freeSequence,
    1327                           update, tableauRow, candidateList);
    1328   }
    1329   const unsigned char *  COIN_RESTRICT internalStatus = model_->internalStatus();
    1330   int *  COIN_RESTRICT index = tableauRow.getIndices();
    1331   double *  COIN_RESTRICT array = tableauRow.denseVector();
    1332   const double * COIN_RESTRICT element = element_;
    1333   const int * COIN_RESTRICT column = column_;
    1334   for (int i=0;i<number;i++) {
    1335     int iRow=piIndex[i];
     1316      update, tableauRow, candidateList);
     1317  }
     1318  const unsigned char *COIN_RESTRICT internalStatus = model_->internalStatus();
     1319  int *COIN_RESTRICT index = tableauRow.getIndices();
     1320  double *COIN_RESTRICT array = tableauRow.denseVector();
     1321  const double *COIN_RESTRICT element = element_;
     1322  const int *COIN_RESTRICT column = column_;
     1323  for (int i = 0; i < number; i++) {
     1324    int iRow = piIndex[i];
    13361325    double piValue = pi[iRow];
    13371326    CoinBigIndex end = rowEnd[iRow];
    1338     for (CoinBigIndex j=rowStart[iRow];j<end;j++) {
     1327    for (CoinBigIndex j = rowStart[iRow]; j < end; j++) {
    13391328      int iColumn = column[j];
    1340       assert (iColumn>=first&&iColumn<last);
     1329      assert(iColumn >= first && iColumn < last);
    13411330      double value = element[j];
    1342       array[iColumn]+= piValue*value;
     1331      array[iColumn] += piValue * value;
    13431332    }
    13441333  }
    13451334  int numberNonZero;
    13461335  int numberRemaining;
    1347   if (iBlock==0) {
     1336  if (iBlock == 0) {
    13481337#if 1
    1349     upperTheta=static_cast<AbcSimplexDual *>(model_)->dualColumn1A();
    1350     numberNonZero=tableauRow.getNumElements(0);
    1351     numberRemaining=candidateList.getNumElements(0);
     1338    upperTheta = static_cast< AbcSimplexDual * >(model_)->dualColumn1A();
     1339    numberNonZero = tableauRow.getNumElements(0);
     1340    numberRemaining = candidateList.getNumElements(0);
    13521341#else
    1353     numberNonZero=0;
    1354     for (int i=0;i<number;i++) {
    1355       int iRow=piIndex[i];
    1356       unsigned char type=internalStatus[iRow];
    1357       if ((type&7)<6) {
    1358         index[numberNonZero]=iRow;
    1359         double piValue=pi[iRow];
    1360         array[numberNonZero++]=piValue;
    1361       }
    1362     }
    1363     numberRemaining=0;
     1342    numberNonZero = 0;
     1343    for (int i = 0; i < number; i++) {
     1344      int iRow = piIndex[i];
     1345      unsigned char type = internalStatus[iRow];
     1346      if ((type & 7) < 6) {
     1347        index[numberNonZero] = iRow;
     1348        double piValue = pi[iRow];
     1349        array[numberNonZero++] = piValue;
     1350      }
     1351    }
     1352    numberRemaining = 0;
    13641353#endif
    13651354  } else {
    1366     numberNonZero=firstIn;
    1367     numberRemaining=firstIn;
    1368   }
    1369   double *  COIN_RESTRICT arrayCandidate=candidateList.denseVector();
    1370   int *  COIN_RESTRICT indexCandidate = candidateList.getIndices();
     1355    numberNonZero = firstIn;
     1356    numberRemaining = firstIn;
     1357  }
     1358  double *COIN_RESTRICT arrayCandidate = candidateList.denseVector();
     1359  int *COIN_RESTRICT indexCandidate = candidateList.getIndices();
    13711360  //printf("first %d last %d firstIn %d lastIn %d\n",
    13721361  //     first,last,firstIn,lastIn);
    1373   const double *  COIN_RESTRICT abcDj = model_->djRegion();
     1362  const double *COIN_RESTRICT abcDj = model_->djRegion();
    13741363  // do first pass to get possibles
    13751364  double bestPossible = 0.0;
     
    13771366  double tentativeTheta = 1.0e25; // try with this much smaller as guess
    13781367  double acceptablePivot = model_->currentAcceptablePivot();
    1379   double dualT=-model_->currentDualTolerance();
    1380   const double multiplier[] = { 1.0, -1.0};
     1368  double dualT = -model_->currentDualTolerance();
     1369  const double multiplier[] = { 1.0, -1.0 };
    13811370  double zeroTolerance = model_->zeroTolerance();
    1382   freeSequence=-1;
     1371  freeSequence = -1;
    13831372  if (model_->ordinaryVariables()) {
    13841373    for (int iSequence = first; iSequence < last; iSequence++) {
    1385       double tableauValue=array[iSequence];
     1374      double tableauValue = array[iSequence];
    13861375      if (tableauValue) {
    1387         array[iSequence]=0.0;
    1388         if (fabs(tableauValue)>zeroTolerance) {
    1389           unsigned char iStatus=internalStatus[iSequence]&7;
    1390           if (iStatus<4) {
    1391             index[numberNonZero]=iSequence;
    1392             array[numberNonZero++]=tableauValue;
    1393             double mult = multiplier[iStatus];
    1394             double alpha = tableauValue * mult;
    1395             double oldValue = abcDj[iSequence] * mult;
    1396             double value = oldValue - tentativeTheta * alpha;
    1397             if (value < dualT) {
    1398               bestPossible = CoinMax(bestPossible, alpha);
    1399               value = oldValue - upperTheta * alpha;
    1400               if (value < dualT && alpha >= acceptablePivot) {
    1401                 upperTheta = (oldValue - dualT) / alpha;
    1402               }
    1403               // add to list
    1404               arrayCandidate[numberRemaining] = alpha;
    1405               indexCandidate[numberRemaining++] = iSequence;
    1406             }
    1407           }
    1408         }
     1376        array[iSequence] = 0.0;
     1377        if (fabs(tableauValue) > zeroTolerance) {
     1378          unsigned char iStatus = internalStatus[iSequence] & 7;
     1379          if (iStatus < 4) {
     1380            index[numberNonZero] = iSequence;
     1381            array[numberNonZero++] = tableauValue;
     1382            double mult = multiplier[iStatus];
     1383            double alpha = tableauValue * mult;
     1384            double oldValue = abcDj[iSequence] * mult;
     1385            double value = oldValue - tentativeTheta * alpha;
     1386            if (value < dualT) {
     1387              bestPossible = CoinMax(bestPossible, alpha);
     1388              value = oldValue - upperTheta * alpha;
     1389              if (value < dualT && alpha >= acceptablePivot) {
     1390                upperTheta = (oldValue - dualT) / alpha;
     1391              }
     1392              // add to list
     1393              arrayCandidate[numberRemaining] = alpha;
     1394              indexCandidate[numberRemaining++] = iSequence;
     1395            }
     1396          }
     1397        }
    14091398      }
    14101399    }
     
    14121401    double badFree = 0.0;
    14131402    double freeAlpha = model_->currentAcceptablePivot();
    1414     int freeSequenceIn=model_->freeSequenceIn();
     1403    int freeSequenceIn = model_->freeSequenceIn();
    14151404    //printf("block %d freeSequence %d acceptable %g\n",iBlock,freeSequenceIn,freeAlpha);
    14161405    double currentDualTolerance = model_->currentDualTolerance();
    14171406    for (int iSequence = first; iSequence < last; iSequence++) {
    1418       double tableauValue=array[iSequence];
     1407      double tableauValue = array[iSequence];
    14191408      if (tableauValue) {
    1420         array[iSequence]=0.0;
    1421         if (fabs(tableauValue)>zeroTolerance) {
    1422           unsigned char iStatus=internalStatus[iSequence]&7;
    1423           if (iStatus<6) {
    1424             if ((iStatus&4)==0) {
    1425               index[numberNonZero]=iSequence;
    1426               array[numberNonZero++]=tableauValue;
    1427               double mult = multiplier[iStatus];
    1428               double alpha = tableauValue * mult;
    1429               double oldValue = abcDj[iSequence] * mult;
    1430               double value = oldValue - tentativeTheta * alpha;
    1431               if (value < dualT) {
    1432                 bestPossible = CoinMax(bestPossible, alpha);
    1433                 value = oldValue - upperTheta * alpha;
    1434                 if (value < dualT && alpha >= acceptablePivot) {
    1435                   upperTheta = (oldValue - dualT) / alpha;
    1436                 }
    1437                 // add to list
    1438                 arrayCandidate[numberRemaining] = alpha;
    1439                 indexCandidate[numberRemaining++] = iSequence;
    1440               }
    1441             } else {
    1442               bool keep;
    1443               index[numberNonZero]=iSequence;
    1444               array[numberNonZero++]=tableauValue;
    1445               bestPossible = CoinMax(bestPossible, fabs(tableauValue));
    1446               double oldValue = abcDj[iSequence];
    1447               // If free has to be very large - should come in via dualRow
    1448               //if (getInternalStatus(iSequence+addSequence)==isFree&&fabs(tableauValue)<1.0e-3)
    1449               //break;
    1450               // may be fake super basic
    1451               if (oldValue > currentDualTolerance) {
    1452                 keep = true;
    1453               } else if (oldValue < -currentDualTolerance) {
    1454                 keep = true;
    1455               } else {
    1456                 if (fabs(tableauValue) > CoinMax(10.0 * acceptablePivot, 1.0e-5)) {
    1457                   keep = true;
    1458                 } else {
    1459                   keep = false;
    1460                   badFree = CoinMax(badFree, fabs(tableauValue));
    1461                 }
    1462               }
     1409        array[iSequence] = 0.0;
     1410        if (fabs(tableauValue) > zeroTolerance) {
     1411          unsigned char iStatus = internalStatus[iSequence] & 7;
     1412          if (iStatus < 6) {
     1413            if ((iStatus & 4) == 0) {
     1414              index[numberNonZero] = iSequence;
     1415              array[numberNonZero++] = tableauValue;
     1416              double mult = multiplier[iStatus];
     1417              double alpha = tableauValue * mult;
     1418              double oldValue = abcDj[iSequence] * mult;
     1419              double value = oldValue - tentativeTheta * alpha;
     1420              if (value < dualT) {
     1421                bestPossible = CoinMax(bestPossible, alpha);
     1422                value = oldValue - upperTheta * alpha;
     1423                if (value < dualT && alpha >= acceptablePivot) {
     1424                  upperTheta = (oldValue - dualT) / alpha;
     1425                }
     1426                // add to list
     1427                arrayCandidate[numberRemaining] = alpha;
     1428                indexCandidate[numberRemaining++] = iSequence;
     1429              }
     1430            } else {
     1431              bool keep;
     1432              index[numberNonZero] = iSequence;
     1433              array[numberNonZero++] = tableauValue;
     1434              bestPossible = CoinMax(bestPossible, fabs(tableauValue));
     1435              double oldValue = abcDj[iSequence];
     1436              // If free has to be very large - should come in via dualRow
     1437              //if (getInternalStatus(iSequence+addSequence)==isFree&&fabs(tableauValue)<1.0e-3)
     1438              //break;
     1439              // may be fake super basic
     1440              if (oldValue > currentDualTolerance) {
     1441                keep = true;
     1442              } else if (oldValue < -currentDualTolerance) {
     1443                keep = true;
     1444              } else {
     1445                if (fabs(tableauValue) > CoinMax(10.0 * acceptablePivot, 1.0e-5)) {
     1446                  keep = true;
     1447                } else {
     1448                  keep = false;
     1449                  badFree = CoinMax(badFree, fabs(tableauValue));
     1450                }
     1451              }
    14631452#if 0
    14641453              if (iSequence==freeSequenceIn)
    14651454                assert (keep);
    14661455#endif
    1467               if (keep) {
     1456              if (keep) {
    14681457#ifdef PAN
    1469                 if (model_->fakeSuperBasic(iSequence)>=0) {
    1470 #endif
    1471                   if (iSequence==freeSequenceIn)
    1472                     tableauValue=COIN_DBL_MAX;
    1473                   // free - choose largest
    1474                   if (fabs(tableauValue) > fabs(freeAlpha)) {
    1475                     freeAlpha = tableauValue;
    1476                     freeSequence = iSequence;
    1477                   }
     1458                if (model_->fakeSuperBasic(iSequence) >= 0) {
     1459#endif
     1460                  if (iSequence == freeSequenceIn)
     1461                    tableauValue = COIN_DBL_MAX;
     1462                  // free - choose largest
     1463                  if (fabs(tableauValue) > fabs(freeAlpha)) {
     1464                    freeAlpha = tableauValue;
     1465                    freeSequence = iSequence;
     1466                  }
    14781467#ifdef PAN
    1479                 }
    1480 #endif
    1481               }
    1482             }
    1483           }
    1484         }
     1468                }
     1469#endif
     1470              }
     1471            }
     1472          }
     1473        }
    14851474      }
    14861475    }
     
    14961485  xxInfo[4][iBlock]=numberRemaining-firstIn;
    14971486#endif
    1498   tableauRow.setNumElementsPartition(iBlock,numberNonZero-firstIn);
    1499   candidateList.setNumElementsPartition(iBlock,numberRemaining-firstIn);
     1487  tableauRow.setNumElementsPartition(iBlock, numberNonZero - firstIn);
     1488  candidateList.setNumElementsPartition(iBlock, numberRemaining - firstIn);
    15001489  return upperTheta;
    15011490}
    15021491// gets sorted tableau row and a possible value of theta
    1503 double 
    1504 AbcMatrix::dualColumn1Row2(double upperTheta, int & freeSequence,
    1505                            const CoinIndexedVector & update,
    1506                            CoinPartitionedVector & tableauRow,
    1507                            CoinPartitionedVector & candidateList) const
     1492double
     1493AbcMatrix::dualColumn1Row2(double upperTheta, int &freeSequence,
     1494  const CoinIndexedVector &update,
     1495  CoinPartitionedVector &tableauRow,
     1496  CoinPartitionedVector &candidateList) const
    15081497{
    15091498  //int first=model_->maximumAbcNumberRows();
    1510   assert(update.getNumElements()==2);
    1511   const double *  COIN_RESTRICT pi=update.denseVector();
    1512   const int *  COIN_RESTRICT piIndex = update.getIndices();
    1513   int *  COIN_RESTRICT index = tableauRow.getIndices();
    1514   double *  COIN_RESTRICT array = tableauRow.denseVector();
    1515   const CoinBigIndex * COIN_RESTRICT rowStart = rowStart_;
    1516   int numberRows=model_->numberRows();
    1517   const CoinBigIndex * COIN_RESTRICT rowEnd = rowStart+numberRows*numberRowBlocks_;
    1518   const double * COIN_RESTRICT element = element_;
    1519   const int * COIN_RESTRICT column = column_;
     1499  assert(update.getNumElements() == 2);
     1500  const double *COIN_RESTRICT pi = update.denseVector();
     1501  const int *COIN_RESTRICT piIndex = update.getIndices();
     1502  int *COIN_RESTRICT index = tableauRow.getIndices();
     1503  double *COIN_RESTRICT array = tableauRow.denseVector();
     1504  const CoinBigIndex *COIN_RESTRICT rowStart = rowStart_;
     1505  int numberRows = model_->numberRows();
     1506  const CoinBigIndex *COIN_RESTRICT rowEnd = rowStart + numberRows * numberRowBlocks_;
     1507  const double *COIN_RESTRICT element = element_;
     1508  const int *COIN_RESTRICT column = column_;
    15201509  int iRow0 = piIndex[0];
    15211510  int iRow1 = piIndex[1];
    15221511  CoinBigIndex end0 = rowEnd[iRow0];
    15231512  CoinBigIndex end1 = rowEnd[iRow1];
    1524   if (end0-rowStart[iRow0]>end1-rowStart[iRow1]) {
    1525     int temp=iRow0;
    1526     iRow0=iRow1;
    1527     iRow1=temp;
     1513  if (end0 - rowStart[iRow0] > end1 - rowStart[iRow1]) {
     1514    int temp = iRow0;
     1515    iRow0 = iRow1;
     1516    iRow1 = temp;
    15281517  }
    15291518  CoinBigIndex start = rowStart[iRow0];
    15301519  CoinBigIndex end = rowEnd[iRow0];
    15311520  double piValue = pi[iRow0];
    1532   double *  COIN_RESTRICT arrayCandidate=candidateList.denseVector();
     1521  double *COIN_RESTRICT arrayCandidate = candidateList.denseVector();
    15331522  int numberNonZero;
    1534   numberNonZero=tableauRow.getNumElements(0);
    1535   int n=numberNonZero;
    1536   for (CoinBigIndex j=start;j<end;j++) {
     1523  numberNonZero = tableauRow.getNumElements(0);
     1524  int n = numberNonZero;
     1525  for (CoinBigIndex j = start; j < end; j++) {
    15371526    int iSequence = column[j];
    15381527    double value = element[j];
    1539     arrayCandidate[iSequence]= piValue*value;
    1540     index[n++]=iSequence;
     1528    arrayCandidate[iSequence] = piValue * value;
     1529    index[n++] = iSequence;
    15411530  }
    15421531  start = rowStart[iRow1];
    15431532  end = rowEnd[iRow1];
    15441533  piValue = pi[iRow1];
    1545   for (CoinBigIndex j=start;j<end;j++) {
     1534  for (CoinBigIndex j = start; j < end; j++) {
    15461535    int iSequence = column[j];
    15471536    double value = element[j];
    15481537    value *= piValue;
    15491538    if (!arrayCandidate[iSequence]) {
    1550       arrayCandidate[iSequence]= value;
    1551       index[n++]=iSequence;
    1552     } else {   
    1553       arrayCandidate[iSequence] += value; 
     1539      arrayCandidate[iSequence] = value;
     1540      index[n++] = iSequence;
     1541    } else {
     1542      arrayCandidate[iSequence] += value;
    15541543    }
    15551544  }
    15561545  // pack down
    1557   const unsigned char *  COIN_RESTRICT internalStatus = model_->internalStatus();
     1546  const unsigned char *COIN_RESTRICT internalStatus = model_->internalStatus();
    15581547  double zeroTolerance = model_->zeroTolerance();
    1559   while (numberNonZero<n) {
    1560     int iSequence=index[numberNonZero];
    1561     double value=arrayCandidate[iSequence];
    1562     arrayCandidate[iSequence]=0.0;
    1563     unsigned char iStatus=internalStatus[iSequence]&7;
    1564     if (fabs(value)>zeroTolerance&&iStatus<6) {
    1565       index[numberNonZero]=iSequence;
    1566       array[numberNonZero++]=value;
     1548  while (numberNonZero < n) {
     1549    int iSequence = index[numberNonZero];
     1550    double value = arrayCandidate[iSequence];
     1551    arrayCandidate[iSequence] = 0.0;
     1552    unsigned char iStatus = internalStatus[iSequence] & 7;
     1553    if (fabs(value) > zeroTolerance && iStatus < 6) {
     1554      index[numberNonZero] = iSequence;
     1555      array[numberNonZero++] = value;
    15671556    } else {
    15681557      // kill
    15691558      n--;
    1570       index[numberNonZero]=index[n];
    1571     }
    1572   }
    1573   tableauRow.setNumElementsPartition(0,numberNonZero);
    1574   return firstPass(model_,0, upperTheta, freeSequence,
    1575                    tableauRow,
    1576                    candidateList) ;
     1559      index[numberNonZero] = index[n];
     1560    }
     1561  }
     1562  tableauRow.setNumElementsPartition(0, numberNonZero);
     1563  return firstPass(model_, 0, upperTheta, freeSequence,
     1564    tableauRow,
     1565    candidateList);
    15771566}
    15781567//static int ixxxxxx=1;
    15791568// gets sorted tableau row and a possible value of theta
    1580 double 
    1581 AbcMatrix::dualColumn1RowFew(int iBlock, double upperTheta, int & freeSequence,
    1582                            const CoinIndexedVector & update,
    1583                            CoinPartitionedVector & tableauRow,
    1584                            CoinPartitionedVector & candidateList) const
     1569double
     1570AbcMatrix::dualColumn1RowFew(int iBlock, double upperTheta, int &freeSequence,
     1571  const CoinIndexedVector &update,
     1572  CoinPartitionedVector &tableauRow,
     1573  CoinPartitionedVector &candidateList) const
    15851574{
    15861575  //int first=model_->maximumAbcNumberRows();
    15871576  int number = update.getNumElements();
    1588   const double *  COIN_RESTRICT pi=update.denseVector();
    1589   const int *  COIN_RESTRICT piIndex = update.getIndices();
    1590   int *  COIN_RESTRICT index = tableauRow.getIndices();
    1591   double *  COIN_RESTRICT array = tableauRow.denseVector();
    1592   const CoinBigIndex * COIN_RESTRICT rowStart = rowStart_;
    1593   int numberRows=model_->numberRows();
    1594   const CoinBigIndex * COIN_RESTRICT rowEnd = rowStart+numberRows*numberRowBlocks_;
    1595   const double * COIN_RESTRICT element = element_;
    1596   const int * COIN_RESTRICT column = column_;
    1597   double *  COIN_RESTRICT arrayCandidate=candidateList.denseVector();
     1577  const double *COIN_RESTRICT pi = update.denseVector();
     1578  const int *COIN_RESTRICT piIndex = update.getIndices();
     1579  int *COIN_RESTRICT index = tableauRow.getIndices();
     1580  double *COIN_RESTRICT array = tableauRow.denseVector();
     1581  const CoinBigIndex *COIN_RESTRICT rowStart = rowStart_;
     1582  int numberRows = model_->numberRows();
     1583  const CoinBigIndex *COIN_RESTRICT rowEnd = rowStart + numberRows * numberRowBlocks_;
     1584  const double *COIN_RESTRICT element = element_;
     1585  const int *COIN_RESTRICT column = column_;
     1586  double *COIN_RESTRICT arrayCandidate = candidateList.denseVector();
    15981587  int numberNonZero;
    1599   assert (iBlock>=0);
    1600   const unsigned char *  COIN_RESTRICT internalStatus = model_->internalStatus();
    1601   int firstIn=blockStart_[iBlock];
    1602   if (iBlock==0) {
    1603     numberNonZero=0;
    1604     for (int i=0;i<number;i++) {
    1605       int iRow=piIndex[i];
    1606       unsigned char type=internalStatus[iRow];
    1607       if ((type&7)<6) {
    1608         index[numberNonZero]=iRow;
    1609         double piValue=pi[iRow];
    1610         array[numberNonZero++]=piValue;
     1588  assert(iBlock >= 0);
     1589  const unsigned char *COIN_RESTRICT internalStatus = model_->internalStatus();
     1590  int firstIn = blockStart_[iBlock];
     1591  if (iBlock == 0) {
     1592    numberNonZero = 0;
     1593    for (int i = 0; i < number; i++) {
     1594      int iRow = piIndex[i];
     1595      unsigned char type = internalStatus[iRow];
     1596      if ((type & 7) < 6) {
     1597        index[numberNonZero] = iRow;
     1598        double piValue = pi[iRow];
     1599        array[numberNonZero++] = piValue;
    16111600      }
    16121601    }
    16131602  } else {
    1614     numberNonZero=firstIn;
    1615   }
    1616   int n=numberNonZero;
    1617   rowStart += iBlock*numberRows;
    1618   rowEnd = rowStart+numberRows;
    1619   for (int i=0;i<number;i++) {
    1620     int iRow=piIndex[i];
     1603    numberNonZero = firstIn;
     1604  }
     1605  int n = numberNonZero;
     1606  rowStart += iBlock * numberRows;
     1607  rowEnd = rowStart + numberRows;
     1608  for (int i = 0; i < number; i++) {
     1609    int iRow = piIndex[i];
    16211610    double piValue = pi[iRow];
    16221611    CoinBigIndex end = rowEnd[iRow];
    1623     for (CoinBigIndex j=rowStart[iRow];j<end;j++) {
     1612    for (CoinBigIndex j = rowStart[iRow]; j < end; j++) {
    16241613      int iColumn = column[j];
    1625       double value = element[j]*piValue;
    1626       double oldValue=arrayCandidate[iColumn];
     1614      double value = element[j] * piValue;
     1615      double oldValue = arrayCandidate[iColumn];
    16271616      value += oldValue;
    16281617      if (!oldValue) {
    1629         arrayCandidate[iColumn]= value;
    1630         index[n++]=iColumn;
    1631       } else if (value) {       
    1632         arrayCandidate[iColumn] = value;
     1618        arrayCandidate[iColumn] = value;
     1619        index[n++] = iColumn;
     1620      } else if (value) {
     1621        arrayCandidate[iColumn] = value;
    16331622      } else {
    1634         arrayCandidate[iColumn] = COIN_INDEXED_REALLY_TINY_ELEMENT;
     1623        arrayCandidate[iColumn] = COIN_INDEXED_REALLY_TINY_ELEMENT;
    16351624      }
    16361625    }
     
    16381627  // pack down
    16391628  double zeroTolerance = model_->zeroTolerance();
    1640   while (numberNonZero<n) {
    1641     int iSequence=index[numberNonZero];
    1642     double value=arrayCandidate[iSequence];
    1643     arrayCandidate[iSequence]=0.0;
    1644     unsigned char iStatus=internalStatus[iSequence]&7;
    1645     if (fabs(value)>zeroTolerance&&iStatus<6) {
    1646       index[numberNonZero]=iSequence;
    1647       array[numberNonZero++]=value;
     1629  while (numberNonZero < n) {
     1630    int iSequence = index[numberNonZero];
     1631    double value = arrayCandidate[iSequence];
     1632    arrayCandidate[iSequence] = 0.0;
     1633    unsigned char iStatus = internalStatus[iSequence] & 7;
     1634    if (fabs(value) > zeroTolerance && iStatus < 6) {
     1635      index[numberNonZero] = iSequence;
     1636      array[numberNonZero++] = value;
    16481637    } else {
    16491638      // kill
    16501639      n--;
    1651       index[numberNonZero]=index[n];
    1652     }
    1653   }
    1654   tableauRow.setNumElementsPartition(iBlock,numberNonZero-firstIn);
    1655   upperTheta = firstPass(model_,iBlock, upperTheta, freeSequence,
    1656                    tableauRow,
    1657                    candidateList) ;
     1640      index[numberNonZero] = index[n];
     1641    }
     1642  }
     1643  tableauRow.setNumElementsPartition(iBlock, numberNonZero - firstIn);
     1644  upperTheta = firstPass(model_, iBlock, upperTheta, freeSequence,
     1645    tableauRow,
     1646    candidateList);
    16581647  return upperTheta;
    16591648}
    16601649// gets sorted tableau row and a possible value of theta
    1661 double 
    1662 AbcMatrix::dualColumn1Row1(double upperTheta, int & freeSequence,
    1663                            const CoinIndexedVector & update,
    1664                            CoinPartitionedVector & tableauRow,
    1665                            CoinPartitionedVector & candidateList) const
    1666 {
    1667   assert(update.getNumElements()==1);
     1650double
     1651AbcMatrix::dualColumn1Row1(double upperTheta, int &freeSequence,
     1652  const CoinIndexedVector &update,
     1653  CoinPartitionedVector &tableauRow,
     1654  CoinPartitionedVector &candidateList) const
     1655{
     1656  assert(update.getNumElements() == 1);
    16681657  int iRow = update.getIndices()[0];
    16691658  double piValue = update.denseVector()[iRow];
    1670   int *  COIN_RESTRICT index = tableauRow.getIndices();
    1671   double *  COIN_RESTRICT array = tableauRow.denseVector();
    1672   const CoinBigIndex * COIN_RESTRICT rowStart = rowStart_;
    1673   int numberRows=model_->numberRows();
    1674   const CoinBigIndex * COIN_RESTRICT rowEnd = rowStart+numberRows*numberRowBlocks_;
     1659  int *COIN_RESTRICT index = tableauRow.getIndices();
     1660  double *COIN_RESTRICT array = tableauRow.denseVector();
     1661  const CoinBigIndex *COIN_RESTRICT rowStart = rowStart_;
     1662  int numberRows = model_->numberRows();
     1663  const CoinBigIndex *COIN_RESTRICT rowEnd = rowStart + numberRows * numberRowBlocks_;
    16751664  CoinBigIndex start = rowStart[iRow];
    16761665  CoinBigIndex end = rowEnd[iRow];
    1677   const double * COIN_RESTRICT element = element_;
    1678   const int * COIN_RESTRICT column = column_;
     1666  const double *COIN_RESTRICT element = element_;
     1667  const int *COIN_RESTRICT column = column_;
    16791668  int numberNonZero;
    16801669  int numberRemaining;
    1681   numberNonZero=tableauRow.getNumElements(0);
     1670  numberNonZero = tableauRow.getNumElements(0);
    16821671  numberRemaining = candidateList.getNumElements(0);
    1683   double *  COIN_RESTRICT arrayCandidate=candidateList.denseVector();
    1684   int *  COIN_RESTRICT indexCandidate = candidateList.getIndices();
    1685   const double *  COIN_RESTRICT abcDj = model_->djRegion();
    1686   const unsigned char *  COIN_RESTRICT internalStatus = model_->internalStatus();
     1672  double *COIN_RESTRICT arrayCandidate = candidateList.denseVector();
     1673  int *COIN_RESTRICT indexCandidate = candidateList.getIndices();
     1674  const double *COIN_RESTRICT abcDj = model_->djRegion();
     1675  const unsigned char *COIN_RESTRICT internalStatus = model_->internalStatus();
    16871676  // do first pass to get possibles
    16881677  double bestPossible = 0.0;
     
    16901679  double tentativeTheta = 1.0e25; // try with this much smaller as guess
    16911680  double acceptablePivot = model_->currentAcceptablePivot();
    1692   double dualT=-model_->currentDualTolerance();
    1693   const double multiplier[] = { 1.0, -1.0};
     1681  double dualT = -model_->currentDualTolerance();
     1682  const double multiplier[] = { 1.0, -1.0 };
    16941683  double zeroTolerance = model_->zeroTolerance();
    1695   freeSequence=-1;
     1684  freeSequence = -1;
    16961685  if (model_->ordinaryVariables()) {
    1697     for (CoinBigIndex j=start;j<end;j++) {
     1686    for (CoinBigIndex j = start; j < end; j++) {
    16981687      int iSequence = column[j];
    16991688      double value = element[j];
    1700       double tableauValue = piValue*value;
    1701       if (fabs(tableauValue)>zeroTolerance) {
    1702         unsigned char iStatus=internalStatus[iSequence]&7;
    1703         if (iStatus<4) {
    1704           index[numberNonZero]=iSequence;
    1705           array[numberNonZero++]=tableauValue;
    1706           double mult = multiplier[iStatus];
    1707           double alpha = tableauValue * mult;
    1708           double oldValue = abcDj[iSequence] * mult;
    1709           double value = oldValue - tentativeTheta * alpha;
    1710           if (value < dualT) {
    1711             bestPossible = CoinMax(bestPossible, alpha);
    1712             value = oldValue - upperTheta * alpha;
    1713             if (value < dualT && alpha >= acceptablePivot) {
    1714               upperTheta = (oldValue - dualT) / alpha;
    1715             }
    1716             // add to list
    1717             arrayCandidate[numberRemaining] = alpha;
    1718             indexCandidate[numberRemaining++] = iSequence;
    1719           }
    1720         }
     1689      double tableauValue = piValue * value;
     1690      if (fabs(tableauValue) > zeroTolerance) {
     1691        unsigned char iStatus = internalStatus[iSequence] & 7;
     1692        if (iStatus < 4) {
     1693          index[numberNonZero] = iSequence;
     1694          array[numberNonZero++] = tableauValue;
     1695          double mult = multiplier[iStatus];
     1696          double alpha = tableauValue * mult;
     1697          double oldValue = abcDj[iSequence] * mult;
     1698          double value = oldValue - tentativeTheta * alpha;
     1699          if (value < dualT) {
     1700            bestPossible = CoinMax(bestPossible, alpha);
     1701            value = oldValue - upperTheta * alpha;
     1702            if (value < dualT && alpha >= acceptablePivot) {
     1703              upperTheta = (oldValue - dualT) / alpha;
     1704            }
     1705            // add to list
     1706            arrayCandidate[numberRemaining] = alpha;
     1707            indexCandidate[numberRemaining++] = iSequence;
     1708          }
     1709        }
    17211710      }
    17221711    }
     
    17241713    double badFree = 0.0;
    17251714    double freeAlpha = model_->currentAcceptablePivot();
    1726     int freeSequenceIn=model_->freeSequenceIn();
     1715    int freeSequenceIn = model_->freeSequenceIn();
    17271716    double currentDualTolerance = model_->currentDualTolerance();
    1728     for (CoinBigIndex j=start;j<end;j++) {
     1717    for (CoinBigIndex j = start; j < end; j++) {
    17291718      int iSequence = column[j];
    17301719      double value = element[j];
    1731       double tableauValue = piValue*value;
    1732       if (fabs(tableauValue)>zeroTolerance) {
    1733         unsigned char iStatus=internalStatus[iSequence]&7;
    1734         if (iStatus<6) {
    1735           if ((iStatus&4)==0) {
    1736             index[numberNonZero]=iSequence;
    1737             array[numberNonZero++]=tableauValue;
    1738             double mult = multiplier[iStatus];
    1739             double alpha = tableauValue * mult;
    1740             double oldValue = abcDj[iSequence] * mult;
    1741             double value = oldValue - tentativeTheta * alpha;
    1742             if (value < dualT) {
    1743               bestPossible = CoinMax(bestPossible, alpha);
    1744               value = oldValue - upperTheta * alpha;
    1745               if (value < dualT && alpha >= acceptablePivot) {
    1746                 upperTheta = (oldValue - dualT) / alpha;
    1747               }
    1748               // add to list
    1749               arrayCandidate[numberRemaining] = alpha;
    1750               indexCandidate[numberRemaining++] = iSequence;
    1751             }
    1752           } else {
    1753             bool keep;
    1754             index[numberNonZero]=iSequence;
    1755             array[numberNonZero++]=tableauValue;
    1756             bestPossible = CoinMax(bestPossible, fabs(tableauValue));
    1757             double oldValue = abcDj[iSequence];
    1758             // If free has to be very large - should come in via dualRow
    1759             //if (getInternalStatus(iSequence+addSequence)==isFree&&fabs(tableauValue)<1.0e-3)
    1760             //break;
    1761             if (oldValue > currentDualTolerance) {
    1762               keep = true;
    1763             } else if (oldValue < -currentDualTolerance) {
    1764               keep = true;
    1765             } else {
    1766               if (fabs(tableauValue) > CoinMax(10.0 * acceptablePivot, 1.0e-5)) {
    1767                 keep = true;
    1768               } else {
    1769                 keep = false;
    1770                 badFree = CoinMax(badFree, fabs(tableauValue));
    1771               }
    1772             }
    1773             if (keep) {
     1720      double tableauValue = piValue * value;
     1721      if (fabs(tableauValue) > zeroTolerance) {
     1722        unsigned char iStatus = internalStatus[iSequence] & 7;
     1723        if (iStatus < 6) {
     1724          if ((iStatus & 4) == 0) {
     1725            index[numberNonZero] = iSequence;
     1726            array[numberNonZero++] = tableauValue;
     1727            double mult = multiplier[iStatus];
     1728            double alpha = tableauValue * mult;
     1729            double oldValue = abcDj[iSequence] * mult;
     1730            double value = oldValue - tentativeTheta * alpha;
     1731            if (value < dualT) {
     1732              bestPossible = CoinMax(bestPossible, alpha);
     1733              value = oldValue - upperTheta * alpha;
     1734              if (value < dualT && alpha >= acceptablePivot) {
     1735                upperTheta = (oldValue - dualT) / alpha;
     1736              }
     1737              // add to list
     1738              arrayCandidate[numberRemaining] = alpha;
     1739              indexCandidate[numberRemaining++] = iSequence;
     1740            }
     1741          } else {
     1742            bool keep;
     1743            index[numberNonZero] = iSequence;
     1744            array[numberNonZero++] = tableauValue;
     1745            bestPossible = CoinMax(bestPossible, fabs(tableauValue));
     1746            double oldValue = abcDj[iSequence];
     1747            // If free has to be very large - should come in via dualRow
     1748            //if (getInternalStatus(iSequence+addSequence)==isFree&&fabs(tableauValue)<1.0e-3)
     1749            //break;
     1750            if (oldValue > currentDualTolerance) {
     1751              keep = true;
     1752            } else if (oldValue < -currentDualTolerance) {
     1753              keep = true;
     1754            } else {
     1755              if (fabs(tableauValue) > CoinMax(10.0 * acceptablePivot, 1.0e-5)) {
     1756                keep = true;
     1757              } else {
     1758                keep = false;
     1759                badFree = CoinMax(badFree, fabs(tableauValue));
     1760              }
     1761            }
     1762            if (keep) {
    17741763#ifdef PAN
    1775               if (model_->fakeSuperBasic(iSequence)>=0) {
    1776 #endif
    1777                 if (iSequence==freeSequenceIn)
    1778                   tableauValue=COIN_DBL_MAX;
    1779                 // free - choose largest
    1780                 if (fabs(tableauValue) > fabs(freeAlpha)) {
    1781                   freeAlpha = tableauValue;
    1782                   freeSequence = iSequence;
    1783                 }
     1764              if (model_->fakeSuperBasic(iSequence) >= 0) {
     1765#endif
     1766                if (iSequence == freeSequenceIn)
     1767                  tableauValue = COIN_DBL_MAX;
     1768                // free - choose largest
     1769                if (fabs(tableauValue) > fabs(freeAlpha)) {
     1770                  freeAlpha = tableauValue;
     1771                  freeSequence = iSequence;
     1772                }
    17841773#ifdef PAN
    1785               }
    1786 #endif
    1787             }
    1788           }
    1789         }
    1790       }
    1791     }
    1792   }
    1793   tableauRow.setNumElementsPartition(0,numberNonZero);
    1794   candidateList.setNumElementsPartition(0,numberRemaining);
     1774              }
     1775#endif
     1776            }
     1777          }
     1778        }
     1779      }
     1780    }
     1781  }
     1782  tableauRow.setNumElementsPartition(0, numberNonZero);
     1783  candidateList.setNumElementsPartition(0, numberRemaining);
    17951784  return upperTheta;
    17961785}
    17971786//#define PARALLEL2
    17981787#ifdef PARALLEL2
    1799 #undef cilk_for 
     1788#undef cilk_for
    18001789#undef cilk_spawn
    18011790#undef cilk_sync
     
    18231812}
    18241813#endif
    1825 void
    1826 AbcMatrix::rebalance() const
     1814void AbcMatrix::rebalance() const
    18271815{
    18281816  int maximumRows = model_->maximumAbcNumberRows();
     
    18361824   */
    18371825#if ABC_PARALLEL
    1838   int howOften=CoinMax(model_->factorization()->maximumPivots(),200);
    1839   if ((model_->numberIterations()%howOften)==0||!startColumnBlock_[1]) {
    1840     int numberCpus=model_->numberCpus();
    1841     if (numberCpus>1) {
    1842       const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts()-maximumRows;
    1843       const unsigned char *  COIN_RESTRICT internalStatus = model_->internalStatus();
    1844       int numberRows=model_->numberRows();
    1845       int total=0;
    1846       for (int iSequence=0;iSequence<numberRows;iSequence++) {
    1847         unsigned char iStatus=internalStatus[iSequence]&7;
    1848         if (iStatus<4)
    1849           total+=3;
    1850       }
    1851       int totalSlacks=total;
     1826  int howOften = CoinMax(model_->factorization()->maximumPivots(), 200);
     1827  if ((model_->numberIterations() % howOften) == 0 || !startColumnBlock_[1]) {
     1828    int numberCpus = model_->numberCpus();
     1829    if (numberCpus > 1) {
     1830      const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts() - maximumRows;
     1831      const unsigned char *COIN_RESTRICT internalStatus = model_->internalStatus();
     1832      int numberRows = model_->numberRows();
     1833      int total = 0;
     1834      for (int iSequence = 0; iSequence < numberRows; iSequence++) {
     1835        unsigned char iStatus = internalStatus[iSequence] & 7;
     1836        if (iStatus < 4)
     1837          total += 3;
     1838      }
     1839      int totalSlacks = total;
    18521840      CoinBigIndex end = columnStart[maximumRows];
    1853       for (int iSequence=maximumRows;iSequence<numberTotal;iSequence++) {
    1854         CoinBigIndex start=end;
    1855         end = columnStart[iSequence+1];
    1856         unsigned char iStatus=internalStatus[iSequence]&7;
    1857         if (iStatus<4)
    1858           total+=5+2*(end-start);
    1859         else
    1860           total+=1;
    1861       }
    1862       int chunk = (total+numberCpus-1)/numberCpus;
    1863       total=totalSlacks;
    1864       int iCpu=0;
    1865       startColumnBlock_[0]=0;
     1841      for (int iSequence = maximumRows; iSequence < numberTotal; iSequence++) {
     1842        CoinBigIndex start = end;
     1843        end = columnStart[iSequence + 1];
     1844        unsigned char iStatus = internalStatus[iSequence] & 7;
     1845        if (iStatus < 4)
     1846          total += 5 + 2 * (end - start);
     1847        else
     1848          total += 1;
     1849      }
     1850      int chunk = (total + numberCpus - 1) / numberCpus;
     1851      total = totalSlacks;
     1852      int iCpu = 0;
     1853      startColumnBlock_[0] = 0;
    18661854      end = columnStart[maximumRows];
    1867       for (int iSequence=maximumRows;iSequence<numberTotal;iSequence++) {
    1868         CoinBigIndex start=end;
    1869         end = columnStart[iSequence+1];
    1870         unsigned char iStatus=internalStatus[iSequence]&7;
    1871         if (iStatus<4)
    1872           total+=5+2*(end-start);
    1873         else
    1874           total+=1;
    1875         if (total>chunk) {
    1876           iCpu++;
    1877           total=0;
    1878           startColumnBlock_[iCpu]=iSequence;
    1879         }
    1880       }
    1881       assert(iCpu<numberCpus);
     1855      for (int iSequence = maximumRows; iSequence < numberTotal; iSequence++) {
     1856        CoinBigIndex start = end;
     1857        end = columnStart[iSequence + 1];
     1858        unsigned char iStatus = internalStatus[iSequence] & 7;
     1859        if (iStatus < 4)
     1860          total += 5 + 2 * (end - start);
     1861        else
     1862          total += 1;
     1863        if (total > chunk) {
     1864          iCpu++;
     1865          total = 0;
     1866          startColumnBlock_[iCpu] = iSequence;
     1867        }
     1868      }
     1869      assert(iCpu < numberCpus);
    18821870      iCpu++;
    1883       for (;iCpu<=numberCpus;iCpu++)
    1884         startColumnBlock_[iCpu]=numberTotal;
    1885       numberColumnBlocks_=numberCpus;
     1871      for (; iCpu <= numberCpus; iCpu++)
     1872        startColumnBlock_[iCpu] = numberTotal;
     1873      numberColumnBlocks_ = numberCpus;
    18861874    } else {
    1887       numberColumnBlocks_=1;
    1888       startColumnBlock_[0]=0;
    1889       startColumnBlock_[1]=numberTotal;
     1875      numberColumnBlocks_ = 1;
     1876      startColumnBlock_[0] = 0;
     1877      startColumnBlock_[1] = numberTotal;
    18901878    }
    18911879  }
    18921880#else
    1893   startColumnBlock_[0]=0;
    1894   startColumnBlock_[1]=numberTotal;
    1895 #endif
    1896 }
    1897 double 
    1898 AbcMatrix::dualColumn1(const CoinIndexedVector & update,
    1899                                   CoinPartitionedVector & tableauRow,CoinPartitionedVector & candidateList) const
     1881  startColumnBlock_[0] = 0;
     1882  startColumnBlock_[1] = numberTotal;
     1883#endif
     1884}
     1885double
     1886AbcMatrix::dualColumn1(const CoinIndexedVector &update,
     1887  CoinPartitionedVector &tableauRow, CoinPartitionedVector &candidateList) const
    19001888{
    19011889  int numberTotal = model_->numberTotal();
     
    19041892  tableauRow.setPackedMode(true);
    19051893  candidateList.setPackedMode(true);
    1906   int number=update.getNumElements();
     1894  int number = update.getNumElements();
    19071895  double upperTheta;
    1908   if (rowStart_&&number<3) {
    1909 #if ABC_INSTRUMENT>1
     1896  if (rowStart_ && number < 3) {
     1897#if ABC_INSTRUMENT > 1
    19101898    {
    1911       int n=update.getNumElements();
    1912       abcPricing[n<19 ? n : 19]++;
     1899      int n = update.getNumElements();
     1900      abcPricing[n < 19 ? n : 19]++;
    19131901    }
    19141902#endif
     
    19161904    // do slacks first
    19171905    int starts[2];
    1918     starts[0]=0;
    1919     starts[1]=numberTotal;
    1920     tableauRow.setPartitions(1,starts);
    1921     candidateList.setPartitions(1,starts);
    1922     upperTheta = static_cast<AbcSimplexDual *>(model_)->dualColumn1A();
     1906    starts[0] = 0;
     1907    starts[1] = numberTotal;
     1908    tableauRow.setPartitions(1, starts);
     1909    candidateList.setPartitions(1, starts);
     1910    upperTheta = static_cast< AbcSimplexDual * >(model_)->dualColumn1A();
    19231911    //assert (upperTheta>-100*model_->dualTolerance()||model_->sequenceIn()>=0||model_->lastFirstFree()>=0);
    1924     int freeSequence=-1;
     1912    int freeSequence = -1;
    19251913    // worth using row copy
    1926     assert (number);
    1927     if (number==2) {
    1928       upperTheta=dualColumn1Row2(upperTheta,freeSequence,update,tableauRow,candidateList);
     1914    assert(number);
     1915    if (number == 2) {
     1916      upperTheta = dualColumn1Row2(upperTheta, freeSequence, update, tableauRow, candidateList);
    19291917    } else {
    1930       upperTheta=dualColumn1Row1(upperTheta,freeSequence,update,tableauRow,candidateList);
    1931     }
    1932     if (freeSequence>=0) {
    1933       int numberNonZero=tableauRow.getNumElements(0);
    1934       const int *  COIN_RESTRICT index = tableauRow.getIndices();
    1935       const double *  COIN_RESTRICT array = tableauRow.denseVector();
     1918      upperTheta = dualColumn1Row1(upperTheta, freeSequence, update, tableauRow, candidateList);
     1919    }
     1920    if (freeSequence >= 0) {
     1921      int numberNonZero = tableauRow.getNumElements(0);
     1922      const int *COIN_RESTRICT index = tableauRow.getIndices();
     1923      const double *COIN_RESTRICT array = tableauRow.denseVector();
    19361924      // search for free coming in
    1937       double freeAlpha=0.0;
    1938       int bestSequence=model_->sequenceIn();
    1939       if (bestSequence>=0)
    1940         freeAlpha=model_->alpha();
     1925      double freeAlpha = 0.0;
     1926      int bestSequence = model_->sequenceIn();
     1927      if (bestSequence >= 0)
     1928        freeAlpha = model_->alpha();
    19411929      index = tableauRow.getIndices();
    19421930      array = tableauRow.denseVector();
    1943       // free variable - search 
    1944       for (int k=0;k<numberNonZero;k++) {
    1945         if (freeSequence==index[k]) {
    1946           if (fabs(freeAlpha)<fabs(array[k])) {
    1947             freeAlpha=array[k];
    1948           }
    1949           break;
    1950         }
    1951       }
    1952       if (model_->sequenceIn()<0||fabs(freeAlpha)>fabs(model_->alpha())) {
    1953         double oldValue = model_->djRegion()[freeSequence];
    1954         model_->setSequenceIn(freeSequence);
    1955         model_->setAlpha(freeAlpha);
    1956         model_->setTheta(oldValue / freeAlpha);
     1931      // free variable - search
     1932      for (int k = 0; k < numberNonZero; k++) {
     1933        if (freeSequence == index[k]) {
     1934          if (fabs(freeAlpha) < fabs(array[k])) {
     1935            freeAlpha = array[k];
     1936          }
     1937          break;
     1938        }
     1939      }
     1940      if (model_->sequenceIn() < 0 || fabs(freeAlpha) > fabs(model_->alpha())) {
     1941        double oldValue = model_->djRegion()[freeSequence];
     1942        model_->setSequenceIn(freeSequence);
     1943        model_->setAlpha(freeAlpha);
     1944        model_->setTheta(oldValue / freeAlpha);
    19571945      }
    19581946    }
     
    19601948    // three or more
    19611949    // need to do better job on dividing up (but wait until vector or by row)
    1962     upperTheta = parallelDual4(static_cast<AbcSimplexDual *>(model_));
     1950    upperTheta = parallelDual4(static_cast< AbcSimplexDual * >(model_));
    19631951  }
    19641952  //tableauRow.compact();
     
    19681956#endif
    19691957  candidateList.computeNumberElements();
    1970   int numberRemaining=candidateList.getNumElements();
    1971   if (!numberRemaining&&model_->sequenceIn()<0) {
     1958  int numberRemaining = candidateList.getNumElements();
     1959  if (!numberRemaining && model_->sequenceIn() < 0) {
    19721960    return COIN_DBL_MAX; // Looks infeasible
    19731961  } else {
     
    19751963  }
    19761964}
    1977 #define _mm256_broadcast_sd(x) static_cast<__m256d> (__builtin_ia32_vbroadcastsd256 (x))
     1965#define _mm256_broadcast_sd(x) static_cast< __m256d >(__builtin_ia32_vbroadcastsd256(x))
    19781966#define _mm256_load_pd(x) *(__m256d *)(x)
    1979 #define _mm256_store_pd (s, x)  *((__m256d *)s)=x
    1980 void
    1981 AbcMatrix::dualColumn1Part(int iBlock, int & sequenceIn, double & upperThetaResult,
    1982                            const CoinIndexedVector & update,
    1983                            CoinPartitionedVector & tableauRow,CoinPartitionedVector & candidateList) const
    1984 {
    1985   double upperTheta=upperThetaResult;
     1967#define _mm256_store_pd (s, x) * ((__m256d *)s) = x
     1968void AbcMatrix::dualColumn1Part(int iBlock, int &sequenceIn, double &upperThetaResult,
     1969  const CoinIndexedVector &update,
     1970  CoinPartitionedVector &tableauRow, CoinPartitionedVector &candidateList) const
     1971{
     1972  double upperTheta = upperThetaResult;
    19861973#if 0
    19871974  double time0=CoinCpuTime();
    19881975#endif
    19891976  int maximumRows = model_->maximumAbcNumberRows();
    1990   int firstIn=startColumnBlock_[iBlock];
    1991   int last = startColumnBlock_[iBlock+1];
     1977  int firstIn = startColumnBlock_[iBlock];
     1978  int last = startColumnBlock_[iBlock + 1];
    19921979  int numberNonZero;
    19931980  int numberRemaining;
    19941981  int first;
    1995   if (firstIn==0) {
    1996     upperTheta=static_cast<AbcSimplexDual *>(model_)->dualColumn1A();
    1997     numberNonZero=tableauRow.getNumElements(0);
     1982  if (firstIn == 0) {
     1983    upperTheta = static_cast< AbcSimplexDual * >(model_)->dualColumn1A();
     1984    numberNonZero = tableauRow.getNumElements(0);
    19981985    numberRemaining = candidateList.getNumElements(0);
    1999     first=maximumRows;
     1986    first = maximumRows;
    20001987  } else {
    2001     first=firstIn;
    2002     numberNonZero=first;
    2003     numberRemaining=first;
    2004   }
    2005   sequenceIn=-1;
     1988    first = firstIn;
     1989    numberNonZero = first;
     1990    numberRemaining = first;
     1991  }
     1992  sequenceIn = -1;
    20061993  // get matrix data pointers
    2007   const int * COIN_RESTRICT row = matrix_->getIndices();
    2008   const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts()-maximumRows;
    2009   const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
    2010   const double *  COIN_RESTRICT pi=update.denseVector();
     1994  const int *COIN_RESTRICT row = matrix_->getIndices();
     1995  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts() - maximumRows;
     1996  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
     1997  const double *COIN_RESTRICT pi = update.denseVector();
    20111998  //const int *  COIN_RESTRICT piIndex = update.getIndices();
    2012   int *  COIN_RESTRICT index = tableauRow.getIndices();
    2013   double *  COIN_RESTRICT array = tableauRow.denseVector();
     1999  int *COIN_RESTRICT index = tableauRow.getIndices();
     2000  double *COIN_RESTRICT array = tableauRow.denseVector();
    20142001  // pivot elements
    2015   double *  COIN_RESTRICT arrayCandidate=candidateList.denseVector();
     2002  double *COIN_RESTRICT arrayCandidate = candidateList.denseVector();
    20162003  // indices
    2017   int *  COIN_RESTRICT indexCandidate = candidateList.getIndices();
    2018   const double *  COIN_RESTRICT abcDj = model_->djRegion();
    2019   const unsigned char *  COIN_RESTRICT internalStatus = model_->internalStatus();
     2004  int *COIN_RESTRICT indexCandidate = candidateList.getIndices();
     2005  const double *COIN_RESTRICT abcDj = model_->djRegion();
     2006  const unsigned char *COIN_RESTRICT internalStatus = model_->internalStatus();
    20202007  // do first pass to get possibles
    20212008  double bestPossible = 0.0;
     
    20232010  double tentativeTheta = 1.0e25; // try with this much smaller as guess
    20242011  double acceptablePivot = model_->currentAcceptablePivot();
    2025   double dualT=-model_->currentDualTolerance();
    2026   const double multiplier[] = { 1.0, -1.0};
     2012  double dualT = -model_->currentDualTolerance();
     2013  const double multiplier[] = { 1.0, -1.0 };
    20272014  double zeroTolerance = model_->zeroTolerance();
    20282015  if (model_->ordinaryVariables()) {
    20292016#ifdef COUNT_COPY
    2030     if (first>maximumRows||last<model_->numberTotal()||false) {
     2017    if (first > maximumRows || last < model_->numberTotal() || false) {
    20312018#endif
    20322019#if 1
    2033     for (int iSequence = first; iSequence < last; iSequence++) {
    2034       unsigned char iStatus=internalStatus[iSequence]&7;
    2035       if (iStatus<4) {
    2036         CoinBigIndex start = columnStart[iSequence];
    2037         CoinBigIndex next = columnStart[iSequence+1];
    2038         double tableauValue = 0.0;
    2039         for (CoinBigIndex j = start; j < next; j++) {
    2040           int jRow = row[j];
    2041           tableauValue += pi[jRow] * elementByColumn[j];
    2042         }
    2043         if (fabs(tableauValue)>zeroTolerance) {
     2020      for (int iSequence = first; iSequence < last; iSequence++) {
     2021        unsigned char iStatus = internalStatus[iSequence] & 7;
     2022        if (iStatus < 4) {
     2023          CoinBigIndex start = columnStart[iSequence];
     2024          CoinBigIndex next = columnStart[iSequence + 1];
     2025          double tableauValue = 0.0;
     2026          for (CoinBigIndex j = start; j < next; j++) {
     2027            int jRow = row[j];
     2028            tableauValue += pi[jRow] * elementByColumn[j];
     2029          }
     2030          if (fabs(tableauValue) > zeroTolerance) {
    20442031#else
    2045     cilk_for (int iSequence = first; iSequence < last; iSequence++) {
    2046       unsigned char iStatus=internalStatus[iSequence]&7;
    2047       if (iStatus<4) {
    2048         CoinBigIndex start = columnStart[iSequence];
    2049         CoinBigIndex next = columnStart[iSequence+1];
    2050         double tableauValue = 0.0;
    2051         for (CoinBigIndex j = start; j < next; j++) {
    2052           int jRow = row[j];
    2053           tableauValue += pi[jRow] * elementByColumn[j];
    2054         }
    2055         array[iSequence]=tableauValue;
     2032    cilk_for(int iSequence = first; iSequence < last; iSequence++)
     2033    {
     2034      unsigned char iStatus = internalStatus[iSequence] & 7;
     2035      if (iStatus < 4) {
     2036        CoinBigIndex start = columnStart[iSequence];
     2037        CoinBigIndex next = columnStart[iSequence + 1];
     2038        double tableauValue = 0.0;
     2039        for (CoinBigIndex j = start; j < next; j++) {
     2040          int jRow = row[j];
     2041          tableauValue += pi[jRow] * elementByColumn[j];
     2042        }
     2043        array[iSequence] = tableauValue;
    20562044      }
    20572045    }
    20582046    //printf("time %.6g\n",CoinCpuTime()-time0);
    20592047    for (int iSequence = first; iSequence < last; iSequence++) {
    2060       double tableauValue=array[iSequence];
     2048      double tableauValue = array[iSequence];
    20612049      if (tableauValue) {
    2062         array[iSequence]=0.0;
    2063         if (fabs(tableauValue)>zeroTolerance) {
    2064           unsigned char iStatus=internalStatus[iSequence]&7;
    2065 #endif
    2066           index[numberNonZero]=iSequence;
    2067           array[numberNonZero++]=tableauValue;
    2068           double mult = multiplier[iStatus];
    2069           double alpha = tableauValue * mult;
    2070           double oldValue = abcDj[iSequence] * mult;
    2071           double value = oldValue - tentativeTheta * alpha;
    2072           if (value < dualT) {
    2073             bestPossible = CoinMax(bestPossible, alpha);
    2074             value = oldValue - upperTheta * alpha;
    2075             if (value < dualT && alpha >= acceptablePivot) {
    2076               upperTheta = (oldValue - dualT) / alpha;
    2077             }
    2078             // add to list
    2079             arrayCandidate[numberRemaining] = alpha;
    2080             indexCandidate[numberRemaining++] = iSequence;
    2081           }
    2082         }
    2083       }
    2084     }
     2050        array[iSequence] = 0.0;
     2051        if (fabs(tableauValue) > zeroTolerance) {
     2052          unsigned char iStatus = internalStatus[iSequence] & 7;
     2053#endif
     2054            index[numberNonZero] = iSequence;
     2055            array[numberNonZero++] = tableauValue;
     2056            double mult = multiplier[iStatus];
     2057            double alpha = tableauValue * mult;
     2058            double oldValue = abcDj[iSequence] * mult;
     2059            double value = oldValue - tentativeTheta * alpha;
     2060            if (value < dualT) {
     2061              bestPossible = CoinMax(bestPossible, alpha);
     2062              value = oldValue - upperTheta * alpha;
     2063              if (value < dualT && alpha >= acceptablePivot) {
     2064                upperTheta = (oldValue - dualT) / alpha;
     2065              }
     2066              // add to list
     2067              arrayCandidate[numberRemaining] = alpha;
     2068              indexCandidate[numberRemaining++] = iSequence;
     2069            }
     2070          }
     2071        }
     2072      }
    20852073#ifdef COUNT_COPY
    20862074    } else {
    2087       const double * COIN_RESTRICT element = countElement_;
    2088       const int * COIN_RESTRICT row = countRow_;
    2089       double temp[4] __attribute__ ((aligned (32)));
    2090       memset(temp,0,sizeof(temp));
    2091       for (int iCount=smallestCount_;iCount<largestCount_;iCount++) {
    2092         int iStart=countFirst_[iCount];
    2093         int number=countFirst_[iCount+1]-iStart;
    2094         if (!number)
    2095           continue;
    2096         const int * COIN_RESTRICT blockRow = row+countStart_[iCount];
    2097         const double * COIN_RESTRICT blockElement = element+countStart_[iCount];
    2098         const int * COIN_RESTRICT sequence = countRealColumn_+iStart;
    2099         // really need to sort and swap to avoid tests
    2100         int numberBlocks=number>>2;
    2101         number &= 3;
    2102         for (int iBlock=0;iBlock<numberBlocks;iBlock++) {
    2103           double tableauValue0=0.0;
    2104           double tableauValue1=0.0;
    2105           double tableauValue2=0.0;
    2106           double tableauValue3=0.0;
    2107           for (int j=0;j<iCount;j++) {
    2108             tableauValue0 += pi[blockRow[0]]*blockElement[0];
    2109             tableauValue1 += pi[blockRow[1]]*blockElement[1];
    2110             tableauValue2 += pi[blockRow[2]]*blockElement[2];
    2111             tableauValue3 += pi[blockRow[3]]*blockElement[3];
    2112             blockRow+=4;
    2113             blockElement+=4;
    2114           }
    2115           if (fabs(tableauValue0)>zeroTolerance) {
    2116             int iSequence=sequence[0];
    2117             unsigned char iStatus=internalStatus[iSequence]&7;
    2118             if (iStatus<4) {
    2119               index[numberNonZero]=iSequence;
    2120               array[numberNonZero++]=tableauValue0;
    2121               double mult = multiplier[iStatus];
    2122               double alpha = tableauValue0 * mult;
    2123               double oldValue = abcDj[iSequence] * mult;
    2124               double value = oldValue - tentativeTheta * alpha;
    2125               if (value < dualT) {
    2126                 bestPossible = CoinMax(bestPossible, alpha);
    2127                 value = oldValue - upperTheta * alpha;
    2128                 if (value < dualT && alpha >= acceptablePivot) {
    2129                   upperTheta = (oldValue - dualT) / alpha;
    2130                 }
    2131                 // add to list
    2132                 arrayCandidate[numberRemaining] = alpha;
    2133                 indexCandidate[numberRemaining++] = iSequence;
    2134               }
    2135             }
    2136           }
    2137           if (fabs(tableauValue1)>zeroTolerance) {
    2138             int iSequence=sequence[1];
    2139             unsigned char iStatus=internalStatus[iSequence]&7;
    2140             if (iStatus<4) {
    2141               index[numberNonZero]=iSequence;
    2142               array[numberNonZero++]=tableauValue1;
    2143               double mult = multiplier[iStatus];
    2144               double alpha = tableauValue1 * mult;
    2145               double oldValue = abcDj[iSequence] * mult;
    2146               double value = oldValue - tentativeTheta * alpha;
    2147               if (value < dualT) {
    2148                 bestPossible = CoinMax(bestPossible, alpha);
    2149                 value = oldValue - upperTheta * alpha;
    2150                 if (value < dualT && alpha >= acceptablePivot) {
    2151                   upperTheta = (oldValue - dualT) / alpha;
    2152                 }
    2153                 // add to list
    2154                 arrayCandidate[numberRemaining] = alpha;
    2155                 indexCandidate[numberRemaining++] = iSequence;
    2156               }
    2157             }
    2158           }
    2159           if (fabs(tableauValue2)>zeroTolerance) {
    2160             int iSequence=sequence[2];
    2161             unsigned char iStatus=internalStatus[iSequence]&7;
    2162             if (iStatus<4) {
    2163               index[numberNonZero]=iSequence;
    2164               array[numberNonZero++]=tableauValue2;
    2165               double mult = multiplier[iStatus];
    2166               double alpha = tableauValue2 * mult;
    2167               double oldValue = abcDj[iSequence] * mult;
    2168               double value = oldValue - tentativeTheta * alpha;
    2169               if (value < dualT) {
    2170                 bestPossible = CoinMax(bestPossible, alpha);
    2171                 value = oldValue - upperTheta * alpha;
    2172                 if (value < dualT && alpha >= acceptablePivot) {
    2173                   upperTheta = (oldValue - dualT) / alpha;
    2174                 }
    2175                 // add to list
    2176                 arrayCandidate[numberRemaining] = alpha;
    2177                 indexCandidate[numberRemaining++] = iSequence;
    2178               }
    2179             }
    2180           }
    2181           if (fabs(tableauValue3)>zeroTolerance) {
    2182             int iSequence=sequence[3];
    2183             unsigned char iStatus=internalStatus[iSequence]&7;
    2184             if (iStatus<4) {
    2185               index[numberNonZero]=iSequence;
    2186               array[numberNonZero++]=tableauValue3;
    2187               double mult = multiplier[iStatus];
    2188               double alpha = tableauValue3 * mult;
    2189               double oldValue = abcDj[iSequence] * mult;
    2190               double value = oldValue - tentativeTheta * alpha;
    2191               if (value < dualT) {
    2192                 bestPossible = CoinMax(bestPossible, alpha);
    2193                 value = oldValue - upperTheta * alpha;
    2194                 if (value < dualT && alpha >= acceptablePivot) {
    2195                   upperTheta = (oldValue - dualT) / alpha;
    2196                 }
    2197                 // add to list
    2198                 arrayCandidate[numberRemaining] = alpha;
    2199                 indexCandidate[numberRemaining++] = iSequence;
    2200               }
    2201             }
    2202           }
    2203           sequence+=4;
    2204         }
    2205         for (int i=0;i<number;i++) {
    2206           int iSequence=sequence[i];
    2207           unsigned char iStatus=internalStatus[iSequence]&7;
    2208           if (iStatus<4) {
    2209             double tableauValue=0.0;
    2210             for (int j=0;j<iCount;j++) {
    2211               int iRow=blockRow[4*j];
    2212               tableauValue += pi[iRow]*blockElement[4*j];
    2213             }
    2214             if (fabs(tableauValue)>zeroTolerance) {
    2215               index[numberNonZero]=iSequence;
    2216               array[numberNonZero++]=tableauValue;
    2217               double mult = multiplier[iStatus];
    2218               double alpha = tableauValue * mult;
    2219               double oldValue = abcDj[iSequence] * mult;
    2220               double value = oldValue - tentativeTheta * alpha;
    2221               if (value < dualT) {
    2222                 bestPossible = CoinMax(bestPossible, alpha);
    2223                 value = oldValue - upperTheta * alpha;
    2224                 if (value < dualT && alpha >= acceptablePivot) {
    2225                   upperTheta = (oldValue - dualT) / alpha;
    2226                 }
    2227                 // add to list
    2228                 arrayCandidate[numberRemaining] = alpha;
    2229                 indexCandidate[numberRemaining++] = iSequence;
    2230               }
    2231             }
    2232           }
    2233           blockRow++;
    2234           blockElement++;
    2235         }
    2236       }
    2237       int numberColumns=model_->numberColumns();
     2075      const double *COIN_RESTRICT element = countElement_;
     2076      const int *COIN_RESTRICT row = countRow_;
     2077      double temp[4] __attribute__((aligned(32)));
     2078      memset(temp, 0, sizeof(temp));
     2079      for (int iCount = smallestCount_; iCount < largestCount_; iCount++) {
     2080        int iStart = countFirst_[iCount];
     2081        int number = countFirst_[iCount + 1] - iStart;
     2082        if (!number)
     2083          continue;
     2084        const int *COIN_RESTRICT blockRow = row + countStart_[iCount];
     2085        const double *COIN_RESTRICT blockElement = element + countStart_[iCount];
     2086        const int *COIN_RESTRICT sequence = countRealColumn_ + iStart;
     2087        // really need to sort and swap to avoid tests
     2088        int numberBlocks = number >> 2;
     2089        number &= 3;
     2090        for (int iBlock = 0; iBlock < numberBlocks; iBlock++) {
     2091          double tableauValue0 = 0.0;
     2092          double tableauValue1 = 0.0;
     2093          double tableauValue2 = 0.0;
     2094          double tableauValue3 = 0.0;
     2095          for (int j = 0; j < iCount; j++) {
     2096            tableauValue0 += pi[blockRow[0]] * blockElement[0];
     2097            tableauValue1 += pi[blockRow[1]] * blockElement[1];
     2098            tableauValue2 += pi[blockRow[2]] * blockElement[2];
     2099            tableauValue3 += pi[blockRow[3]] * blockElement[3];
     2100            blockRow += 4;
     2101            blockElement += 4;
     2102          }
     2103          if (fabs(tableauValue0) > zeroTolerance) {
     2104            int iSequence = sequence[0];
     2105            unsigned char iStatus = internalStatus[iSequence] & 7;
     2106            if (iStatus < 4) {
     2107              index[numberNonZero] = iSequence;
     2108              array[numberNonZero++] = tableauValue0;
     2109              double mult = multiplier[iStatus];
     2110              double alpha = tableauValue0 * mult;
     2111              double oldValue = abcDj[iSequence] * mult;
     2112              double value = oldValue - tentativeTheta * alpha;
     2113              if (value < dualT) {
     2114                bestPossible = CoinMax(bestPossible, alpha);
     2115                value = oldValue - upperTheta * alpha;
     2116                if (value < dualT && alpha >= acceptablePivot) {
     2117                  upperTheta = (oldValue - dualT) / alpha;
     2118                }
     2119                // add to list
     2120                arrayCandidate[numberRemaining] = alpha;
     2121                indexCandidate[numberRemaining++] = iSequence;
     2122              }
     2123            }
     2124          }
     2125          if (fabs(tableauValue1) > zeroTolerance) {
     2126            int iSequence = sequence[1];
     2127            unsigned char iStatus = internalStatus[iSequence] & 7;
     2128            if (iStatus < 4) {
     2129              index[numberNonZero] = iSequence;
     2130              array[numberNonZero++] = tableauValue1;
     2131              double mult = multiplier[iStatus];
     2132              double alpha = tableauValue1 * mult;
     2133              double oldValue = abcDj[iSequence] * mult;
     2134              double value = oldValue - tentativeTheta * alpha;
     2135              if (value < dualT) {
     2136                bestPossible = CoinMax(bestPossible, alpha);
     2137                value = oldValue - upperTheta * alpha;
     2138                if (value < dualT && alpha >= acceptablePivot) {
     2139                  upperTheta = (oldValue - dualT) / alpha;
     2140                }
     2141                // add to list
     2142                arrayCandidate[numberRemaining] = alpha;
     2143                indexCandidate[numberRemaining++] = iSequence;
     2144              }
     2145            }
     2146          }
     2147          if (fabs(tableauValue2) > zeroTolerance) {
     2148            int iSequence = sequence[2];
     2149            unsigned char iStatus = internalStatus[iSequence] & 7;
     2150            if (iStatus < 4) {
     2151              index[numberNonZero] = iSequence;
     2152              array[numberNonZero++] = tableauValue2;
     2153              double mult = multiplier[iStatus];
     2154              double alpha = tableauValue2 * mult;
     2155              double oldValue = abcDj[iSequence] * mult;
     2156              double value = oldValue - tentativeTheta * alpha;
     2157              if (value < dualT) {
     2158                bestPossible = CoinMax(bestPossible, alpha);
     2159                value = oldValue - upperTheta * alpha;
     2160                if (value < dualT && alpha >= acceptablePivot) {
     2161                  upperTheta = (oldValue - dualT) / alpha;
     2162                }
     2163                // add to list
     2164                arrayCandidate[numberRemaining] = alpha;
     2165                indexCandidate[numberRemaining++] = iSequence;
     2166              }
     2167            }
     2168          }
     2169          if (fabs(tableauValue3) > zeroTolerance) {
     2170            int iSequence = sequence[3];
     2171            unsigned char iStatus = internalStatus[iSequence] & 7;
     2172            if (iStatus < 4) {
     2173              index[numberNonZero] = iSequence;
     2174              array[numberNonZero++] = tableauValue3;
     2175              double mult = multiplier[iStatus];
     2176              double alpha = tableauValue3 * mult;
     2177              double oldValue = abcDj[iSequence] * mult;
     2178              double value = oldValue - tentativeTheta * alpha;
     2179              if (value < dualT) {
     2180                bestPossible = CoinMax(bestPossible, alpha);
     2181                value = oldValue - upperTheta * alpha;
     2182                if (value < dualT && alpha >= acceptablePivot) {
     2183                  upperTheta = (oldValue - dualT) / alpha;
     2184                }
     2185                // add to list
     2186                arrayCandidate[numberRemaining] = alpha;
     2187                indexCandidate[numberRemaining++] = iSequence;
     2188              }
     2189            }
     2190          }
     2191          sequence += 4;
     2192        }
     2193        for (int i = 0; i < number; i++) {
     2194          int iSequence = sequence[i];
     2195          unsigned char iStatus = internalStatus[iSequence] & 7;
     2196          if (iStatus < 4) {
     2197            double tableauValue = 0.0;
     2198            for (int j = 0; j < iCount; j++) {
     2199              int iRow = blockRow[4 * j];
     2200              tableauValue += pi[iRow] * blockElement[4 * j];
     2201            }
     2202            if (fabs(tableauValue) > zeroTolerance) {
     2203              index[numberNonZero] = iSequence;
     2204              array[numberNonZero++] = tableauValue;
     2205              double mult = multiplier[iStatus];
     2206              double alpha = tableauValue * mult;
     2207              double oldValue = abcDj[iSequence] * mult;
     2208              double value = oldValue - tentativeTheta * alpha;
     2209              if (value < dualT) {
     2210                bestPossible = CoinMax(bestPossible, alpha);
     2211                value = oldValue - upperTheta * alpha;
     2212                if (value < dualT && alpha >= acceptablePivot) {
     2213                  upperTheta = (oldValue - dualT) / alpha;
     2214                }
     2215                // add to list
     2216                arrayCandidate[numberRemaining] = alpha;
     2217                indexCandidate[numberRemaining++] = iSequence;
     2218              }
     2219            }
     2220          }
     2221          blockRow++;
     2222          blockElement++;
     2223        }
     2224      }
     2225      int numberColumns = model_->numberColumns();
    22382226      // large ones
    2239       const CoinBigIndex * COIN_RESTRICT largeStart = countStartLarge_-countFirst_[MAX_COUNT];
    2240       for (int i=countFirst_[MAX_COUNT];i<numberColumns;i++) {
    2241         int iSequence=countRealColumn_[i];
    2242         unsigned char iStatus=internalStatus[iSequence]&7;
    2243         if (iStatus<4) {
    2244           CoinBigIndex start = largeStart[i];
    2245           CoinBigIndex next = largeStart[i+1];
    2246           double tableauValue = 0.0;
    2247           for (CoinBigIndex j = start; j < next; j++) {
    2248             int jRow = row[j];
    2249             tableauValue += pi[jRow] * element[j];
    2250           }
    2251           if (fabs(tableauValue)>zeroTolerance) {
    2252             index[numberNonZero]=iSequence;
    2253             array[numberNonZero++]=tableauValue;
    2254             double mult = multiplier[iStatus];
    2255             double alpha = tableauValue * mult;
    2256             double oldValue = abcDj[iSequence] * mult;
    2257             double value = oldValue - tentativeTheta * alpha;
    2258             if (value < dualT) {
    2259               bestPossible = CoinMax(bestPossible, alpha);
    2260               value = oldValue - upperTheta * alpha;
    2261               if (value < dualT && alpha >= acceptablePivot) {
    2262                 upperTheta = (oldValue - dualT) / alpha;
    2263               }
    2264               // add to list
    2265               arrayCandidate[numberRemaining] = alpha;
    2266               indexCandidate[numberRemaining++] = iSequence;
    2267             }
    2268           }
    2269         }
     2227      const CoinBigIndex *COIN_RESTRICT largeStart = countStartLarge_ - countFirst_[MAX_COUNT];
     2228      for (int i = countFirst_[MAX_COUNT]; i < numberColumns; i++) {
     2229        int iSequence = countRealColumn_[i];
     2230        unsigned char iStatus = internalStatus[iSequence] & 7;
     2231        if (iStatus < 4) {
     2232          CoinBigIndex start = largeStart[i];
     2233          CoinBigIndex next = largeStart[i + 1];
     2234          double tableauValue = 0.0;
     2235          for (CoinBigIndex j = start; j < next; j++) {
     2236            int jRow = row[j];
     2237            tableauValue += pi[jRow] * element[j];
     2238          }
     2239          if (fabs(tableauValue) > zeroTolerance) {
     2240            index[numberNonZero] = iSequence;
     2241            array[numberNonZero++] = tableauValue;
     2242            double mult = multiplier[iStatus];
     2243            double alpha = tableauValue * mult;
     2244            double oldValue = abcDj[iSequence] * mult;
     2245            double value = oldValue - tentativeTheta * alpha;
     2246            if (value < dualT) {
     2247              bestPossible = CoinMax(bestPossible, alpha);
     2248              value = oldValue - upperTheta * alpha;
     2249              if (value < dualT && alpha >= acceptablePivot) {
     2250                upperTheta = (oldValue - dualT) / alpha;
     2251              }
     2252              // add to list
     2253              arrayCandidate[numberRemaining] = alpha;
     2254              indexCandidate[numberRemaining++] = iSequence;
     2255            }
     2256          }
     2257        }
    22702258      }
    22712259    }
     
    22742262    double badFree = 0.0;
    22752263    double freePivot = model_->currentAcceptablePivot();
    2276     int freeSequenceIn=model_->freeSequenceIn();
     2264    int freeSequenceIn = model_->freeSequenceIn();
    22772265    double currentDualTolerance = model_->currentDualTolerance();
    22782266    for (int iSequence = first; iSequence < last; iSequence++) {
    2279       unsigned char iStatus=internalStatus[iSequence]&7;
    2280       if (iStatus<6) {
    2281         CoinBigIndex start = columnStart[iSequence];
    2282         CoinBigIndex next = columnStart[iSequence+1];
    2283         double tableauValue = 0.0;
    2284         for (CoinBigIndex j = start; j < next; j++) {
    2285           int jRow = row[j];
    2286           tableauValue += pi[jRow] * elementByColumn[j];
    2287         }
    2288         if (fabs(tableauValue)>zeroTolerance) {
    2289           if ((iStatus&4)==0) {
    2290             index[numberNonZero]=iSequence;
    2291             array[numberNonZero++]=tableauValue;
    2292             double mult = multiplier[iStatus];
    2293             double alpha = tableauValue * mult;
    2294             double oldValue = abcDj[iSequence] * mult;
    2295             double value = oldValue - tentativeTheta * alpha;
    2296             if (value < dualT) {
    2297               bestPossible = CoinMax(bestPossible, alpha);
    2298               value = oldValue - upperTheta * alpha;
    2299               if (value < dualT && alpha >= acceptablePivot) {
    2300                 upperTheta = (oldValue - dualT) / alpha;
    2301               }
    2302               // add to list
    2303               arrayCandidate[numberRemaining] = alpha;
    2304               indexCandidate[numberRemaining++] = iSequence;
    2305             }
    2306           } else {
    2307             bool keep;
    2308             index[numberNonZero]=iSequence;
    2309             array[numberNonZero++]=tableauValue;
    2310             bestPossible = CoinMax(bestPossible, fabs(tableauValue));
    2311             double oldValue = abcDj[iSequence];
    2312             // If free has to be very large - should come in via dualRow
    2313             //if (getInternalStatus(iSequence+addSequence)==isFree&&fabs(tableauValue)<1.0e-3)
    2314             //break;
    2315             if (oldValue > currentDualTolerance) {
    2316               keep = true;
    2317             } else if (oldValue < -currentDualTolerance) {
    2318               keep = true;
    2319             } else {
    2320               if (fabs(tableauValue) > CoinMax(10.0 * acceptablePivot, 1.0e-5)) {
    2321                 keep = true;
    2322               } else {
    2323                 keep = false;
    2324                 badFree = CoinMax(badFree, fabs(tableauValue));
    2325               }
    2326             }
    2327             if (keep) {
     2267      unsigned char iStatus = internalStatus[iSequence] & 7;
     2268      if (iStatus < 6) {
     2269        CoinBigIndex start = columnStart[iSequence];
     2270        CoinBigIndex next = columnStart[iSequence + 1];
     2271        double tableauValue = 0.0;
     2272        for (CoinBigIndex j = start; j < next; j++) {
     2273          int jRow = row[j];
     2274          tableauValue += pi[jRow] * elementByColumn[j];
     2275        }
     2276        if (fabs(tableauValue) > zeroTolerance) {
     2277          if ((iStatus & 4) == 0) {
     2278            index[numberNonZero] = iSequence;
     2279            array[numberNonZero++] = tableauValue;
     2280            double mult = multiplier[iStatus];
     2281            double alpha = tableauValue * mult;
     2282            double oldValue = abcDj[iSequence] * mult;
     2283            double value = oldValue - tentativeTheta * alpha;
     2284            if (value < dualT) {
     2285              bestPossible = CoinMax(bestPossible, alpha);
     2286              value = oldValue - upperTheta * alpha;
     2287              if (value < dualT && alpha >= acceptablePivot) {
     2288                upperTheta = (oldValue - dualT) / alpha;
     2289              }
     2290              // add to list
     2291              arrayCandidate[numberRemaining] = alpha;
     2292              indexCandidate[numberRemaining++] = iSequence;
     2293            }
     2294          } else {
     2295            bool keep;
     2296            index[numberNonZero] = iSequence;
     2297            array[numberNonZero++] = tableauValue;
     2298            bestPossible = CoinMax(bestPossible, fabs(tableauValue));
     2299            double oldValue = abcDj[iSequence];
     2300            // If free has to be very large - should come in via dualRow
     2301            //if (getInternalStatus(iSequence+addSequence)==isFree&&fabs(tableauValue)<1.0e-3)
     2302            //break;
     2303            if (oldValue > currentDualTolerance) {
     2304              keep = true;
     2305            } else if (oldValue < -currentDualTolerance) {
     2306              keep = true;
     2307            } else {
     2308              if (fabs(tableauValue) > CoinMax(10.0 * acceptablePivot, 1.0e-5)) {
     2309                keep = true;
     2310              } else {
     2311                keep = false;
     2312                badFree = CoinMax(badFree, fabs(tableauValue));
     2313              }
     2314            }
     2315            if (keep) {
    23282316#ifdef PAN
    2329               if (model_->fakeSuperBasic(iSequence)>=0) {
    2330 #endif
    2331                 if (iSequence==freeSequenceIn)
    2332                   tableauValue=COIN_DBL_MAX;
    2333                 // free - choose largest
    2334                 if (fabs(tableauValue) > fabs(freePivot)) {
    2335                   freePivot = tableauValue;
    2336                   sequenceIn = iSequence;
    2337                 }
     2317              if (model_->fakeSuperBasic(iSequence) >= 0) {
     2318#endif
     2319                if (iSequence == freeSequenceIn)
     2320                  tableauValue = COIN_DBL_MAX;
     2321                // free - choose largest
     2322                if (fabs(tableauValue) > fabs(freePivot)) {
     2323                  freePivot = tableauValue;
     2324                  sequenceIn = iSequence;
     2325                }
    23382326#ifdef PAN
    2339               }
    2340 #endif
    2341             }
    2342           }
    2343         }
    2344       }
    2345     }
    2346   } 
     2327              }
     2328#endif
     2329            }
     2330          }
     2331        }
     2332      }
     2333    }
     2334  }
    23472335  // adjust
    23482336  numberNonZero -= firstIn;
    23492337  numberRemaining -= firstIn;
    2350   tableauRow.setNumElementsPartition(iBlock,numberNonZero);
    2351   candidateList.setNumElementsPartition(iBlock,numberRemaining);
    2352   if (!numberRemaining&&model_->sequenceIn()<0) {
    2353  
    2354     upperThetaResult=COIN_DBL_MAX; // Looks infeasible
     2338  tableauRow.setNumElementsPartition(iBlock, numberNonZero);
     2339  candidateList.setNumElementsPartition(iBlock, numberRemaining);
     2340  if (!numberRemaining && model_->sequenceIn() < 0) {
     2341
     2342    upperThetaResult = COIN_DBL_MAX; // Looks infeasible
    23552343  } else {
    2356     upperThetaResult=upperTheta;
    2357   }
    2358 }
    2359 // gets tableau row - returns number of slacks in block
    2360 int
    2361 AbcMatrix::primalColumnRow(int iBlock,bool doByRow,const CoinIndexedVector & updates,
    2362                            CoinPartitionedVector & tableauRow) const
     2344    upperThetaResult = upperTheta;
     2345  }
     2346}
     2347// gets tableau row - returns number of slacks in block
     2348int AbcMatrix::primalColumnRow(int iBlock, bool doByRow, const CoinIndexedVector &updates,
     2349  CoinPartitionedVector &tableauRow) const
    23632350{
    23642351  int maximumRows = model_->maximumAbcNumberRows();
    2365   int first=tableauRow.startPartition(iBlock);
    2366   int last=tableauRow.startPartition(iBlock+1);
    2367   if (doByRow) { 
    2368     assert(first==blockStart_[iBlock]);
    2369     assert(last==blockStart_[iBlock+1]);
     2352  int first = tableauRow.startPartition(iBlock);
     2353  int last = tableauRow.startPartition(iBlock + 1);
     2354  if (doByRow) {
     2355    assert(first == blockStart_[iBlock]);
     2356    assert(last == blockStart_[iBlock + 1]);
    23702357  } else {
    2371     assert(first==startColumnBlock_[iBlock]);
    2372     assert(last==startColumnBlock_[iBlock+1]);
    2373   }
    2374   int numberSlacks=updates.getNumElements();
    2375   int numberNonZero=0;
    2376   const double *  COIN_RESTRICT pi=updates.denseVector();
    2377   const int *  COIN_RESTRICT piIndex = updates.getIndices();
    2378   int *  COIN_RESTRICT index = tableauRow.getIndices()+first;
    2379   double *  COIN_RESTRICT array = tableauRow.denseVector()+first;
    2380   const unsigned char *  COIN_RESTRICT internalStatus = model_->internalStatus();
     2358    assert(first == startColumnBlock_[iBlock]);
     2359    assert(last == startColumnBlock_[iBlock + 1]);
     2360  }
     2361  int numberSlacks = updates.getNumElements();
     2362  int numberNonZero = 0;
     2363  const double *COIN_RESTRICT pi = updates.denseVector();
     2364  const int *COIN_RESTRICT piIndex = updates.getIndices();
     2365  int *COIN_RESTRICT index = tableauRow.getIndices() + first;
     2366  double *COIN_RESTRICT array = tableauRow.denseVector() + first;
     2367  const unsigned char *COIN_RESTRICT internalStatus = model_->internalStatus();
    23812368  if (!iBlock) {
    2382     numberNonZero=numberSlacks;
    2383     for (int i=0;i<numberSlacks;i++) {
    2384       int iRow=piIndex[i];
    2385       index[i]=iRow;
    2386       array[i]=pi[iRow]; // ? what if small or basic
    2387     }
    2388     first=maximumRows;
     2369    numberNonZero = numberSlacks;
     2370    for (int i = 0; i < numberSlacks; i++) {
     2371      int iRow = piIndex[i];
     2372      index[i] = iRow;
     2373      array[i] = pi[iRow]; // ? what if small or basic
     2374    }
     2375    first = maximumRows;
    23892376  }
    23902377  double zeroTolerance = model_->zeroTolerance();
    23912378  if (doByRow) {
    2392     int numberRows=model_->numberRows();
    2393     const CoinBigIndex * COIN_RESTRICT rowStart = rowStart_+iBlock*numberRows;
    2394     const CoinBigIndex * COIN_RESTRICT rowEnd = rowStart+numberRows;
    2395     const double * COIN_RESTRICT element = element_;
    2396     const int * COIN_RESTRICT column = column_;
    2397     if (numberSlacks>1) {
    2398       double *  COIN_RESTRICT arrayTemp = tableauRow.denseVector();
    2399       for (int i=0;i<numberSlacks;i++) {
    2400         int iRow=piIndex[i];
    2401         double piValue = pi[iRow];
    2402         CoinBigIndex end = rowEnd[iRow];
    2403         for (CoinBigIndex j=rowStart[iRow];j<end;j++) {
    2404           int iColumn = column[j];
    2405           arrayTemp[iColumn] += element[j]*piValue;
    2406         }
    2407       }
    2408       for (int iSequence=first;iSequence<last;iSequence++) {
    2409         double tableauValue = arrayTemp[iSequence];
    2410         if (tableauValue) {
    2411           arrayTemp[iSequence]=0.0;
    2412           if (fabs(tableauValue)>zeroTolerance) {
    2413             unsigned char iStatus=internalStatus[iSequence]&7;
    2414             if (iStatus<6) {
    2415               index[numberNonZero]=iSequence;
    2416               array[numberNonZero++]=tableauValue;
    2417             }
    2418           }
    2419         }
     2379    int numberRows = model_->numberRows();
     2380    const CoinBigIndex *COIN_RESTRICT rowStart = rowStart_ + iBlock * numberRows;
     2381    const CoinBigIndex *COIN_RESTRICT rowEnd = rowStart + numberRows;
     2382    const double *COIN_RESTRICT element = element_;
     2383    const int *COIN_RESTRICT column = column_;
     2384    if (numberSlacks > 1) {
     2385      double *COIN_RESTRICT arrayTemp = tableauRow.denseVector();
     2386      for (int i = 0; i < numberSlacks; i++) {
     2387        int iRow = piIndex[i];
     2388        double piValue = pi[iRow];
     2389        CoinBigIndex end = rowEnd[iRow];
     2390        for (CoinBigIndex j = rowStart[iRow]; j < end; j++) {
     2391          int iColumn = column[j];
     2392          arrayTemp[iColumn] += element[j] * piValue;
     2393        }
     2394      }
     2395      for (int iSequence = first; iSequence < last; iSequence++) {
     2396        double tableauValue = arrayTemp[iSequence];
     2397        if (tableauValue) {
     2398          arrayTemp[iSequence] = 0.0;
     2399          if (fabs(tableauValue) > zeroTolerance) {
     2400            unsigned char iStatus = internalStatus[iSequence] & 7;
     2401            if (iStatus < 6) {
     2402              index[numberNonZero] = iSequence;
     2403              array[numberNonZero++] = tableauValue;
     2404            }
     2405          }
     2406        }
    24202407      }
    24212408    } else {
    2422       int iRow=piIndex[0];
     2409      int iRow = piIndex[0];
    24232410      double piValue = pi[iRow];
    24242411      CoinBigIndex end = rowEnd[iRow];
    2425       for (CoinBigIndex j=rowStart[iRow];j<end;j++) {
    2426         int iSequence = column[j];
    2427         double tableauValue = element[j]*piValue;
    2428         if (fabs(tableauValue)>zeroTolerance) {
    2429           unsigned char iStatus=internalStatus[iSequence]&7;
    2430           if (iStatus<6) {
    2431             index[numberNonZero]=iSequence;
    2432             array[numberNonZero++]=tableauValue;
    2433           }
    2434         }
     2412      for (CoinBigIndex j = rowStart[iRow]; j < end; j++) {
     2413        int iSequence = column[j];
     2414        double tableauValue = element[j] * piValue;
     2415        if (fabs(tableauValue) > zeroTolerance) {
     2416          unsigned char iStatus = internalStatus[iSequence] & 7;
     2417          if (iStatus < 6) {
     2418            index[numberNonZero] = iSequence;
     2419            array[numberNonZero++] = tableauValue;
     2420          }
     2421        }
    24352422      }
    24362423    }
    24372424  } else {
    24382425    // get matrix data pointers
    2439     const int * COIN_RESTRICT row = matrix_->getIndices();
    2440     const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts()-maximumRows;
    2441     const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
    2442     const double *  COIN_RESTRICT pi=updates.denseVector();
     2426    const int *COIN_RESTRICT row = matrix_->getIndices();
     2427    const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts() - maximumRows;
     2428    const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
     2429    const double *COIN_RESTRICT pi = updates.denseVector();
    24432430    for (int iSequence = first; iSequence < last; iSequence++) {
    2444       unsigned char iStatus=internalStatus[iSequence]&7;
    2445       if (iStatus<6) {
    2446         CoinBigIndex start = columnStart[iSequence];
    2447         CoinBigIndex next = columnStart[iSequence+1];
    2448         double tableauValue = 0.0;
    2449         for (CoinBigIndex j = start; j < next; j++) {
    2450           int jRow = row[j];
    2451           tableauValue += pi[jRow] * elementByColumn[j];
    2452         }
    2453         if (fabs(tableauValue)>zeroTolerance) {
    2454           index[numberNonZero]=iSequence;
    2455           array[numberNonZero++]=tableauValue;
    2456         }
    2457       }
    2458     }
    2459   }
    2460   tableauRow.setNumElementsPartition(iBlock,numberNonZero);
     2431      unsigned char iStatus = internalStatus[iSequence] & 7;
     2432      if (iStatus < 6) {
     2433        CoinBigIndex start = columnStart[iSequence];
     2434        CoinBigIndex next = columnStart[iSequence + 1];
     2435        double tableauValue = 0.0;
     2436        for (CoinBigIndex j = start; j < next; j++) {
     2437          int jRow = row[j];
     2438          tableauValue += pi[jRow] * elementByColumn[j];
     2439        }
     2440        if (fabs(tableauValue) > zeroTolerance) {
     2441          index[numberNonZero] = iSequence;
     2442          array[numberNonZero++] = tableauValue;
     2443        }
     2444      }
     2445    }
     2446  }
     2447  tableauRow.setNumElementsPartition(iBlock, numberNonZero);
    24612448  return numberSlacks;
    24622449}
    24632450// Get sequenceIn when Dantzig
    2464  int
    2465 AbcMatrix::pivotColumnDantzig(const CoinIndexedVector & updates,
    2466                               CoinPartitionedVector & spare) const
    2467  {
    2468    int maximumRows = model_->maximumAbcNumberRows();
    2469    // rebalance
    2470    rebalance();
    2471    spare.setPackedMode(true);
    2472    bool useRowCopy=false;
    2473    if (rowStart_) {
    2474      // see if worth using row copy
    2475      int number=updates.getNumElements();
    2476      assert (number);
    2477      useRowCopy = (number<<2)<maximumRows;
    2478    }
    2479    const int * starts;
    2480    if (useRowCopy)
    2481      starts=blockStart();
    2482    else
    2483      starts=startColumnBlock();
     2451int AbcMatrix::pivotColumnDantzig(const CoinIndexedVector &updates,
     2452  CoinPartitionedVector &spare) const
     2453{
     2454  int maximumRows = model_->maximumAbcNumberRows();
     2455  // rebalance
     2456  rebalance();
     2457  spare.setPackedMode(true);
     2458  bool useRowCopy = false;
     2459  if (rowStart_) {
     2460    // see if worth using row copy
     2461    int number = updates.getNumElements();
     2462    assert(number);
     2463    useRowCopy = (number << 2) < maximumRows;
     2464  }
     2465  const int *starts;
     2466  if (useRowCopy)
     2467    starts = blockStart();
     2468  else
     2469    starts = startColumnBlock();
    24842470#if ABC_PARALLEL
    24852471#define NUMBER_BLOCKS NUMBER_ROW_BLOCKS
    2486   int numberBlocks=CoinMin(NUMBER_BLOCKS,model_->numberCpus());
     2472  int numberBlocks = CoinMin(NUMBER_BLOCKS, model_->numberCpus());
    24872473#else
    24882474#define NUMBER_BLOCKS 1
    2489   int numberBlocks=NUMBER_BLOCKS;
     2475  int numberBlocks = NUMBER_BLOCKS;
    24902476#endif
    24912477#if ABC_PARALLEL
    2492   if (model_->parallelMode()==0)
    2493     numberBlocks=1;
    2494 #endif
    2495    spare.setPartitions(numberBlocks,starts);
    2496    int which[NUMBER_BLOCKS];
    2497    double best[NUMBER_BLOCKS];
    2498    for (int i=0;i<numberBlocks-1;i++)
    2499      which[i]=cilk_spawn pivotColumnDantzig(i,useRowCopy,updates,spare,best[i]);
    2500    which[numberBlocks-1]=pivotColumnDantzig(numberBlocks-1,useRowCopy,updates,
    2501                                                        spare,best[numberBlocks-1]);
    2502    cilk_sync;
    2503    int bestSequence=-1;
    2504    double bestValue=model_->dualTolerance();
    2505    for (int i=0;i<numberBlocks;i++) {
    2506      if (best[i]>bestValue) {
    2507        bestValue=best[i];
    2508        bestSequence=which[i];
    2509      }
    2510    }
    2511    return bestSequence;
    2512 }
    2513  // Get sequenceIn when Dantzig (One block)
    2514  int
    2515    AbcMatrix::pivotColumnDantzig(int iBlock, bool doByRow,const CoinIndexedVector & updates,
    2516                                  CoinPartitionedVector & spare,
    2517                               double & bestValue) const
    2518  {
    2519    // could rewrite to do more directly
    2520    int numberSlacks=primalColumnRow(iBlock,doByRow,updates,spare);
    2521    double * COIN_RESTRICT reducedCost = model_->djRegion();
    2522    int first=spare.startPartition(iBlock);
    2523    int last=spare.startPartition(iBlock+1);
    2524    int *  COIN_RESTRICT index = spare.getIndices()+first;
    2525    double *  COIN_RESTRICT array = spare.denseVector()+first;
    2526    int numberNonZero=spare.getNumElements(iBlock);
    2527    int bestSequence=-1;
    2528    double bestDj=model_->dualTolerance();
    2529    double bestFreeDj = model_->dualTolerance();
    2530    int bestFreeSequence = -1;
    2531    // redo LOTS so signs for slacks and columns same
    2532    if (!first) {
    2533      first = model_->maximumAbcNumberRows();
    2534      for (int i=0;i<numberSlacks;i++) {
    2535        int iSequence = index[i];
    2536        double value = reducedCost[iSequence];
    2537        value += array[i];
    2538        array[i] = 0.0;
    2539        reducedCost[iSequence] = value;
    2540      }
    2541 #ifndef CLP_PRIMAL_SLACK_MULTIPLIER
     2478  if (model_->parallelMode() == 0)
     2479    numberBlocks = 1;
     2480#endif
     2481  spare.setPartitions(numberBlocks, starts);
     2482  int which[NUMBER_BLOCKS];
     2483  double best[NUMBER_BLOCKS];
     2484  for (int i = 0; i < numberBlocks - 1; i++)
     2485    which[i] = cilk_spawn pivotColumnDantzig(i, useRowCopy, updates, spare, best[i]);
     2486  which[numberBlocks - 1] = pivotColumnDantzig(numberBlocks - 1, useRowCopy, updates,
     2487    spare, best[numberBlocks - 1]);
     2488  cilk_sync;
     2489  int bestSequence = -1;
     2490  double bestValue = model_->dualTolerance();
     2491  for (int i = 0; i < numberBlocks; i++) {
     2492    if (best[i] > bestValue) {
     2493      bestValue = best[i];
     2494      bestSequence = which[i];
     2495    }
     2496  }
     2497  return bestSequence;
     2498}
     2499// Get sequenceIn when Dantzig (One block)
     2500int AbcMatrix::pivotColumnDantzig(int iBlock, bool doByRow, const CoinIndexedVector &updates,
     2501  CoinPartitionedVector &spare,
     2502  double &bestValue) const
     2503{
     2504  // could rewrite to do more directly
     2505  int numberSlacks = primalColumnRow(iBlock, doByRow, updates, spare);
     2506  double *COIN_RESTRICT reducedCost = model_->djRegion();
     2507  int first = spare.startPartition(iBlock);
     2508  int last = spare.startPartition(iBlock + 1);
     2509  int *COIN_RESTRICT index = spare.getIndices() + first;
     2510  double *COIN_RESTRICT array = spare.denseVector() + first;
     2511  int numberNonZero = spare.getNumElements(iBlock);
     2512  int bestSequence = -1;
     2513  double bestDj = model_->dualTolerance();
     2514  double bestFreeDj = model_->dualTolerance();
     2515  int bestFreeSequence = -1;
     2516  // redo LOTS so signs for slacks and columns same
     2517  if (!first) {
     2518    first = model_->maximumAbcNumberRows();
     2519    for (int i = 0; i < numberSlacks; i++) {
     2520      int iSequence = index[i];
     2521      double value = reducedCost[iSequence];
     2522      value += array[i];
     2523      array[i] = 0.0;
     2524      reducedCost[iSequence] = value;
     2525    }
     2526#ifndef CLP_PRIMAL_SLACK_MULTIPLIER
    25422527#define CLP_PRIMAL_SLACK_MULTIPLIER 1.0
    25432528#endif
    2544      int numberRows=model_->numberRows();
    2545      for (int iSequence=0 ; iSequence < numberRows; iSequence++) {
    2546        // check flagged variable
    2547        if (!model_->flagged(iSequence)) {
    2548          double value = reducedCost[iSequence] * CLP_PRIMAL_SLACK_MULTIPLIER;
    2549          AbcSimplex::Status status = model_->getInternalStatus(iSequence);
    2550          
    2551          switch(status) {
    2552            
    2553          case AbcSimplex::basic:
    2554          case AbcSimplex::isFixed:
    2555            break;
    2556          case AbcSimplex::isFree:
    2557          case AbcSimplex::superBasic:
    2558            if (fabs(value) > bestFreeDj) {
    2559              bestFreeDj = fabs(value);
    2560              bestFreeSequence = iSequence;
    2561            }
    2562            break;
    2563          case AbcSimplex::atUpperBound:
    2564            if (value > bestDj) {
    2565              bestDj = value;
    2566              bestSequence = iSequence;
    2567            }
    2568            break;
    2569          case AbcSimplex::atLowerBound:
    2570            if (value < -bestDj) {
    2571              bestDj = -value;
    2572              bestSequence = iSequence;
    2573            }
    2574          }
    2575        }
    2576      }
    2577    }
    2578    for (int i=numberSlacks;i<numberNonZero;i++) {
    2579      int iSequence = index[i];
    2580      double value = reducedCost[iSequence];
    2581      value += array[i];
    2582      array[i] = 0.0;
    2583      reducedCost[iSequence] = value;
    2584    }
    2585    for (int iSequence=first ; iSequence < last; iSequence++) {
    2586      // check flagged variable
    2587      if (!model_->flagged(iSequence)) {
    2588        double value = reducedCost[iSequence];
    2589        AbcSimplex::Status status = model_->getInternalStatus(iSequence);
    2590        
    2591        switch(status) {
    2592          
    2593        case AbcSimplex::basic:
    2594        case AbcSimplex::isFixed:
    2595          break;
    2596        case AbcSimplex::isFree:
    2597        case AbcSimplex::superBasic:
    2598          if (fabs(value) > bestFreeDj) {
    2599            bestFreeDj = fabs(value);
    2600            bestFreeSequence = iSequence;
    2601          }
    2602          break;
    2603        case AbcSimplex::atUpperBound:
    2604          if (value > bestDj) {
    2605            bestDj = value;
    2606            bestSequence = iSequence;
    2607          }
    2608          break;
    2609        case AbcSimplex::atLowerBound:
    2610          if (value < -bestDj) {
    2611            bestDj = -value;
    2612            bestSequence = iSequence;
    2613          }
    2614        }
    2615      }
    2616    }
    2617    spare.setNumElementsPartition(iBlock,0);
    2618    // bias towards free
    2619    if (bestFreeSequence >= 0 && bestFreeDj > 0.1 * bestDj) {
    2620      bestSequence = bestFreeSequence;
    2621      bestDj=10.0*bestFreeDj;
    2622    }
    2623    bestValue=bestDj;
    2624    return bestSequence;
    2625  }
    2626 // gets tableau row and dj row - returns number of slacks in block
    2627  int
    2628    AbcMatrix::primalColumnRowAndDjs(int iBlock,const CoinIndexedVector & updateTableau,
    2629                                     const CoinIndexedVector & updateDjs,
    2630                                     CoinPartitionedVector & tableauRow) const
     2529    int numberRows = model_->numberRows();
     2530    for (int iSequence = 0; iSequence < numberRows; iSequence++) {
     2531      // check flagged variable
     2532      if (!model_->flagged(iSequence)) {
     2533        double value = reducedCost[iSequence] * CLP_PRIMAL_SLACK_MULTIPLIER;
     2534        AbcSimplex::Status status = model_->getInternalStatus(iSequence);
     2535
     2536        switch (status) {
     2537
     2538        case AbcSimplex::basic:
     2539        case AbcSimplex::isFixed:
     2540          break;
     2541        case AbcSimplex::isFree:
     2542        case AbcSimplex::superBasic:
     2543          if (fabs(value) > bestFreeDj) {
     2544            bestFreeDj = fabs(value);
     2545            bestFreeSequence = iSequence;
     2546          }
     2547          break;
     2548        case AbcSimplex::atUpperBound:
     2549          if (value > bestDj) {
     2550            bestDj = value;
     2551            bestSequence = iSequence;
     2552          }
     2553          break;
     2554        case AbcSimplex::atLowerBound:
     2555          if (value < -bestDj) {
     2556            bestDj = -value;
     2557            bestSequence = iSequence;
     2558          }
     2559        }
     2560      }
     2561    }
     2562  }
     2563  for (int i = numberSlacks; i < numberNonZero; i++) {
     2564    int iSequence = index[i];
     2565    double value = reducedCost[iSequence];
     2566    value += array[i];
     2567    array[i] = 0.0;
     2568    reducedCost[iSequence] = value;
     2569  }
     2570  for (int iSequence = first; iSequence < last; iSequence++) {
     2571    // check flagged variable
     2572    if (!model_->flagged(iSequence)) {
     2573      double value = reducedCost[iSequence];
     2574      AbcSimplex::Status status = model_->getInternalStatus(iSequence);
     2575
     2576      switch (status) {
     2577
     2578      case AbcSimplex::basic:
     2579      case AbcSimplex::isFixed:
     2580        break;
     2581      case AbcSimplex::isFree:
     2582      case AbcSimplex::superBasic:
     2583        if (fabs(value) > bestFreeDj) {
     2584          bestFreeDj = fabs(value);
     2585          bestFreeSequence = iSequence;
     2586        }
     2587        break;
     2588      case AbcSimplex::atUpperBound:
     2589        if (value > bestDj) {
     2590          bestDj = value;
     2591          bestSequence = iSequence;
     2592        }
     2593        break;
     2594      case AbcSimplex::atLowerBound:
     2595        if (value < -bestDj) {
     2596          bestDj = -value;
     2597          bestSequence = iSequence;
     2598        }
     2599      }
     2600    }
     2601  }
     2602  spare.setNumElementsPartition(iBlock, 0);
     2603  // bias towards free
     2604  if (bestFreeSequence >= 0 && bestFreeDj > 0.1 * bestDj) {
     2605    bestSequence = bestFreeSequence;
     2606    bestDj = 10.0 * bestFreeDj;
     2607  }
     2608  bestValue = bestDj;
     2609  return bestSequence;
     2610}
     2611// gets tableau row and dj row - returns number of slacks in block
     2612int AbcMatrix::primalColumnRowAndDjs(int iBlock, const CoinIndexedVector &updateTableau,
     2613  const CoinIndexedVector &updateDjs,
     2614  CoinPartitionedVector &tableauRow) const
    26312615{
    26322616  int maximumRows = model_->maximumAbcNumberRows();
    2633   int first=tableauRow.startPartition(iBlock);
    2634   int last=tableauRow.startPartition(iBlock+1);
    2635   assert(first==startColumnBlock_[iBlock]);
    2636   assert(last==startColumnBlock_[iBlock+1]);
    2637   int numberSlacks=updateTableau.getNumElements();
    2638   int numberNonZero=0;
    2639   const double *  COIN_RESTRICT piTableau=updateTableau.denseVector();
    2640   const double *  COIN_RESTRICT pi=updateDjs.denseVector();
    2641   int *  COIN_RESTRICT index = tableauRow.getIndices()+first;
    2642   double *  COIN_RESTRICT array = tableauRow.denseVector()+first;
    2643   const unsigned char *  COIN_RESTRICT internalStatus = model_->internalStatus();
    2644   double * COIN_RESTRICT reducedCost = model_->djRegion();
     2617  int first = tableauRow.startPartition(iBlock);
     2618  int last = tableauRow.startPartition(iBlock + 1);
     2619  assert(first == startColumnBlock_[iBlock]);
     2620  assert(last == startColumnBlock_[iBlock + 1]);
     2621  int numberSlacks = updateTableau.getNumElements();
     2622  int numberNonZero = 0;
     2623  const double *COIN_RESTRICT piTableau = updateTableau.denseVector();
     2624  const double *COIN_RESTRICT pi = updateDjs.denseVector();
     2625  int *COIN_RESTRICT index = tableauRow.getIndices() + first;
     2626  double *COIN_RESTRICT array = tableauRow.denseVector() + first;
     2627  const unsigned char *COIN_RESTRICT internalStatus = model_->internalStatus();
     2628  double *COIN_RESTRICT reducedCost = model_->djRegion();
    26452629  if (!iBlock) {
    2646     const int *  COIN_RESTRICT piIndexTableau = updateTableau.getIndices();
    2647     numberNonZero=numberSlacks;
    2648     for (int i=0;i<numberSlacks;i++) {
    2649       int iRow=piIndexTableau[i];
    2650       index[i]=iRow;
    2651       array[i]=piTableau[iRow]; // ? what if small or basic
    2652     }
    2653     const int *  COIN_RESTRICT piIndex = updateDjs.getIndices();
    2654     int number=updateDjs.getNumElements();
    2655     for (int i=0;i<number;i++) {
    2656       int iRow=piIndex[i];
     2630    const int *COIN_RESTRICT piIndexTableau = updateTableau.getIndices();
     2631    numberNonZero = numberSlacks;
     2632    for (int i = 0; i < numberSlacks; i++) {
     2633      int iRow = piIndexTableau[i];
     2634      index[i] = iRow;
     2635      array[i] = piTableau[iRow]; // ? what if small or basic
     2636    }
     2637    const int *COIN_RESTRICT piIndex = updateDjs.getIndices();
     2638    int number = updateDjs.getNumElements();
     2639    for (int i = 0; i < number; i++) {
     2640      int iRow = piIndex[i];
    26572641      reducedCost[iRow] -= pi[iRow]; // sign?
    26582642    }
    2659     first=maximumRows;
     2643    first = maximumRows;
    26602644  }
    26612645  double zeroTolerance = model_->zeroTolerance();
    26622646  // get matrix data pointers
    2663   const int * COIN_RESTRICT row = matrix_->getIndices();
    2664   const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts()-maximumRows;
    2665   const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
     2647  const int *COIN_RESTRICT row = matrix_->getIndices();
     2648  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts() - maximumRows;
     2649  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
    26662650  for (int iSequence = first; iSequence < last; iSequence++) {
    2667     unsigned char iStatus=internalStatus[iSequence]&7;
    2668     if (iStatus<6) {
     2651    unsigned char iStatus = internalStatus[iSequence] & 7;
     2652    if (iStatus < 6) {
    26692653      CoinBigIndex start = columnStart[iSequence];
    2670       CoinBigIndex next = columnStart[iSequence+1];
     2654      CoinBigIndex next = columnStart[iSequence + 1];
    26712655      double tableauValue = 0.0;
    2672       double djValue=reducedCost[iSequence];
     2656      double djValue = reducedCost[iSequence];
    26732657      for (CoinBigIndex j = start; j < next; j++) {
    2674         int jRow = row[j];
    2675         tableauValue += piTableau[jRow] * elementByColumn[j];
    2676         djValue -= pi[jRow] * elementByColumn[j]; // sign?
    2677       }
    2678       reducedCost[iSequence]=djValue;
    2679       if (fabs(tableauValue)>zeroTolerance) {
    2680         index[numberNonZero]=iSequence;
    2681         array[numberNonZero++]=tableauValue;
    2682       }
    2683     }
    2684   }
    2685   tableauRow.setNumElementsPartition(iBlock,numberNonZero);
     2658        int jRow = row[j];
     2659        tableauValue += piTableau[jRow] * elementByColumn[j];
     2660        djValue -= pi[jRow] * elementByColumn[j]; // sign?
     2661      }
     2662      reducedCost[iSequence] = djValue;
     2663      if (fabs(tableauValue) > zeroTolerance) {
     2664        index[numberNonZero] = iSequence;
     2665        array[numberNonZero++] = tableauValue;
     2666      }
     2667    }
     2668  }
     2669  tableauRow.setNumElementsPartition(iBlock, numberNonZero);
    26862670  return numberSlacks;
    26872671}
     
    26902674   otherwise use updateForDjs
    26912675*/
    2692 int
    2693 AbcMatrix::primalColumnDouble(int iBlock,CoinPartitionedVector & updateForTableauRow,
    2694                               CoinPartitionedVector & updateForDjs,
    2695                               const CoinIndexedVector & updateForWeights,
    2696                               CoinPartitionedVector & spareColumn1,
    2697                               double * infeasibilities,
    2698                               double referenceIn, double devex,
    2699                               // Array for exact devex to say what is in reference framework
    2700                               unsigned int * reference,
    2701                               double * weights, double scaleFactor) const
     2676int AbcMatrix::primalColumnDouble(int iBlock, CoinPartitionedVector &updateForTableauRow,
     2677  CoinPartitionedVector &updateForDjs,
     2678  const CoinIndexedVector &updateForWeights,
     2679  CoinPartitionedVector &spareColumn1,
     2680  double *infeasibilities,
     2681  double referenceIn, double devex,
     2682  // Array for exact devex to say what is in reference framework
     2683  unsigned int *reference,
     2684  double *weights, double scaleFactor) const
    27022685{
    27032686  int maximumRows = model_->maximumAbcNumberRows();
    2704   int first=startColumnBlock_[iBlock];
    2705   int last=startColumnBlock_[iBlock+1];
    2706   const double *  COIN_RESTRICT piTableau=updateForTableauRow.denseVector();
    2707   const double *  COIN_RESTRICT pi=updateForDjs.denseVector();
    2708   const double *  COIN_RESTRICT piWeights=updateForWeights.denseVector();
    2709   const unsigned char *  COIN_RESTRICT internalStatus = model_->internalStatus();
    2710   double * COIN_RESTRICT reducedCost = model_->djRegion();
     2687  int first = startColumnBlock_[iBlock];
     2688  int last = startColumnBlock_[iBlock + 1];
     2689  const double *COIN_RESTRICT piTableau = updateForTableauRow.denseVector();
     2690  const double *COIN_RESTRICT pi = updateForDjs.denseVector();
     2691  const double *COIN_RESTRICT piWeights = updateForWeights.denseVector();
     2692  const unsigned char *COIN_RESTRICT internalStatus = model_->internalStatus();
     2693  double *COIN_RESTRICT reducedCost = model_->djRegion();
    27112694  double tolerance = model_->currentDualTolerance();
    27122695  // use spareColumn to track new infeasibilities
    2713   int * COIN_RESTRICT newInfeas = spareColumn1.getIndices()+first;
    2714   int numberNewInfeas=0;
     2696  int *COIN_RESTRICT newInfeas = spareColumn1.getIndices() + first;
     2697  int numberNewInfeas = 0;
    27152698  // we can't really trust infeasibilities if there is dual error
    27162699  // this coding has to mimic coding in checkDualSolution
    27172700  double error = CoinMin(1.0e-2, model_->largestDualError());
    27182701  // allow tolerance at least slightly bigger than standard
    2719   tolerance = tolerance +  error;
    2720   double mult[2]={1.0,-1.0};
    2721   bool updateWeights=devex!=0.0;
    2722 #define isReference(i) (((reference[i>>5] >> (i & 31)) & 1) != 0)
     2702  tolerance = tolerance + error;
     2703  double mult[2] = { 1.0, -1.0 };
     2704  bool updateWeights = devex != 0.0;
     2705#define isReference(i) (((reference[i >> 5] >> (i & 31)) & 1) != 0)
    27232706  // Scale factor is other way round so tableau row is 1.0* while dj update is scaleFactor*
    27242707  if (!iBlock) {
    2725     int numberSlacks=updateForTableauRow.getNumElements();
    2726     const int *  COIN_RESTRICT piIndexTableau = updateForTableauRow.getIndices();
     2708    int numberSlacks = updateForTableauRow.getNumElements();
     2709    const int *COIN_RESTRICT piIndexTableau = updateForTableauRow.getIndices();
    27272710    if (!scaleFactor) {
    2728       const int *  COIN_RESTRICT piIndex = updateForDjs.getIndices();
    2729       int number=updateForDjs.getNumElements();
    2730       for (int i=0;i<number;i++) {
    2731         int iRow=piIndex[i];
    2732         int iStatus=internalStatus[iRow]&7;
    2733         double value=reducedCost[iRow];
    2734         value += pi[iRow];
    2735         if (iStatus<4) {
    2736           reducedCost[iRow]=value;
    2737           value *= mult[iStatus];
    2738           if (value<0.0) {
    2739             if (!infeasibilities[iRow])
    2740               newInfeas[numberNewInfeas++]=iRow;
    2741 #ifdef CLP_PRIMAL_SLACK_MULTIPLIER 
    2742             infeasibilities[iRow]=value*value*CLP_PRIMAL_SLACK_MULTIPLIER;
     2711      const int *COIN_RESTRICT piIndex = updateForDjs.getIndices();
     2712      int number = updateForDjs.getNumElements();
     2713      for (int i = 0; i < number; i++) {
     2714        int iRow = piIndex[i];
     2715        int iStatus = internalStatus[iRow] & 7;
     2716        double value = reducedCost[iRow];
     2717        value += pi[iRow];
     2718        if (iStatus < 4) {
     2719          reducedCost[iRow] = value;
     2720          value *= mult[iStatus];
     2721          if (value < 0.0) {
     2722            if (!infeasibilities[iRow])
     2723              newInfeas[numberNewInfeas++] = iRow;
     2724#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
     2725            infeasibilities[iRow] = value * value * CLP_PRIMAL_SLACK_MULTIPLIER;
    27432726#else
    2744             infeasibilities[iRow]=value*value;
    2745 #endif
    2746           } else {
    2747             if (infeasibilities[iRow])
    2748               infeasibilities[iRow]=COIN_INDEXED_REALLY_TINY_ELEMENT;
    2749           }
    2750         }
    2751       }
    2752     } else if (scaleFactor!=COIN_DBL_MAX) {
    2753       for (int i=0;i<numberSlacks;i++) {
    2754         int iRow=piIndexTableau[i];
    2755         int iStatus=internalStatus[iRow]&7;
    2756         if (iStatus<4) {
    2757           double value=reducedCost[iRow];
    2758           value += scaleFactor*piTableau[iRow]; // sign?
    2759           reducedCost[iRow]=value;
    2760           value *= mult[iStatus];
    2761           if (value<0.0) {
    2762             if (!infeasibilities[iRow])
    2763               newInfeas[numberNewInfeas++]=iRow;
    2764 #ifdef CLP_PRIMAL_SLACK_MULTIPLIER 
    2765             infeasibilities[iRow]=value*value*CLP_PRIMAL_SLACK_MULTIPLIER;
     2727            infeasibilities[iRow] = value * value;
     2728#endif
     2729          } else {
     2730            if (infeasibilities[iRow])
     2731              infeasibilities[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
     2732          }
     2733        }
     2734      }
     2735    } else if (scaleFactor != COIN_DBL_MAX) {
     2736      for (int i = 0; i < numberSlacks; i++) {
     2737        int iRow = piIndexTableau[i];
     2738        int iStatus = internalStatus[iRow] & 7;
     2739        if (iStatus < 4) {
     2740          double value = reducedCost[iRow];
     2741          value += scaleFactor * piTableau[iRow]; // sign?
     2742          reducedCost[iRow] = value;
     2743          value *= mult[iStatus];
     2744          if (value < 0.0) {
     2745            if (!infeasibilities[iRow])
     2746              newInfeas[numberNewInfeas++] = iRow;
     2747#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
     2748            infeasibilities[iRow] = value * value * CLP_PRIMAL_SLACK_MULTIPLIER;
    27662749#else
    2767             infeasibilities[iRow]=value*value;
    2768 #endif
    2769           } else {
    2770             if (infeasibilities[iRow])
    2771               infeasibilities[iRow]=COIN_INDEXED_REALLY_TINY_ELEMENT;
    2772           }
    2773         }
     2750            infeasibilities[iRow] = value * value;
     2751#endif
     2752          } else {
     2753            if (infeasibilities[iRow])
     2754              infeasibilities[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
     2755          }
     2756        }
    27742757      }
    27752758    }
    27762759    if (updateWeights) {
    2777       for (int i=0;i<numberSlacks;i++) {
    2778         int iRow=piIndexTableau[i];
    2779         double modification=piWeights[iRow];
    2780         double thisWeight = weights[iRow];
    2781         double pivot = piTableau[iRow];
    2782         double pivotSquared = pivot * pivot;
    2783         thisWeight += pivotSquared * devex - pivot * modification;
    2784         if (thisWeight < DEVEX_TRY_NORM) {
    2785           if (referenceIn < 0.0) {
    2786             // steepest
    2787             thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
    2788           } else {
    2789             // exact
    2790             thisWeight = referenceIn * pivotSquared;
    2791             if (isReference(iRow))
    2792               thisWeight += 1.0;
    2793             thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
    2794           }
    2795         }
    2796         weights[iRow] = thisWeight;
    2797       }
    2798     }
    2799     first=maximumRows;
     2760      for (int i = 0; i < numberSlacks; i++) {
     2761        int iRow = piIndexTableau[i];
     2762        double modification = piWeights[iRow];
     2763        double thisWeight = weights[iRow];
     2764        double pivot = piTableau[iRow];
     2765        double pivotSquared = pivot * pivot;
     2766        thisWeight += pivotSquared * devex - pivot * modification;
     2767        if (thisWeight < DEVEX_TRY_NORM) {
     2768          if (referenceIn < 0.0) {
     2769            // steepest
     2770            thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
     2771          } else {
     2772            // exact
     2773            thisWeight = referenceIn * pivotSquared;
     2774            if (isReference(iRow))
     2775              thisWeight += 1.0;
     2776            thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
     2777          }
     2778        }
     2779        weights[iRow] = thisWeight;
     2780      }
     2781    }
     2782    first = maximumRows;
    28002783  }
    28012784  double zeroTolerance = model_->zeroTolerance();
    2802   double freeTolerance = FREE_ACCEPT * tolerance;;
     2785  double freeTolerance = FREE_ACCEPT * tolerance;
     2786  ;
    28032787  // get matrix data pointers
    2804   const int * COIN_RESTRICT row = matrix_->getIndices();
    2805   const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts()-maximumRows;
    2806   const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
    2807   bool updateDjs=scaleFactor!=COIN_DBL_MAX;
     2788  const int *COIN_RESTRICT row = matrix_->getIndices();
     2789  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts() - maximumRows;
     2790  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
     2791  bool updateDjs = scaleFactor != COIN_DBL_MAX;
    28082792  for (int iSequence = first; iSequence < last; iSequence++) {
    2809     unsigned char iStatus=internalStatus[iSequence]&7;
    2810     if (iStatus<6) {
     2793    unsigned char iStatus = internalStatus[iSequence] & 7;
     2794    if (iStatus < 6) {
    28112795      CoinBigIndex start = columnStart[iSequence];
    2812       CoinBigIndex next = columnStart[iSequence+1];
     2796      CoinBigIndex next = columnStart[iSequence + 1];
    28132797      double tableauValue = 0.0;
    2814       double djValue=reducedCost[iSequence];
     2798      double djValue = reducedCost[iSequence];
    28152799      if (!scaleFactor) {
    2816         for (CoinBigIndex j = start; j < next; j++) {
    2817           int jRow = row[j];
    2818           tableauValue += piTableau[jRow] * elementByColumn[j];
    2819           djValue += pi[jRow] * elementByColumn[j];
    2820         }
     2800        for (CoinBigIndex j = start; j < next; j++) {
     2801          int jRow = row[j];
     2802          tableauValue += piTableau[jRow] * elementByColumn[j];
     2803          djValue += pi[jRow] * elementByColumn[j];
     2804        }
    28212805      } else {
    2822         for (CoinBigIndex j = start; j < next; j++) {
    2823           int jRow = row[j];
    2824           tableauValue += piTableau[jRow] * elementByColumn[j];
    2825         }
    2826         djValue += tableauValue*scaleFactor; // sign?
     2806        for (CoinBigIndex j = start; j < next; j++) {
     2807          int jRow = row[j];
     2808          tableauValue += piTableau[jRow] * elementByColumn[j];
     2809        }
     2810        djValue += tableauValue * scaleFactor; // sign?
    28272811      }
    28282812      if (updateDjs) {
    2829         reducedCost[iSequence]=djValue;
    2830         if (iStatus<4) {
    2831           djValue *= mult[iStatus];
    2832           if (djValue<0.0) {
    2833             if (!infeasibilities[iSequence])
    2834               newInfeas[numberNewInfeas++]=iSequence;
    2835             infeasibilities[iSequence]=djValue*djValue;
    2836           } else {
    2837             if (infeasibilities[iSequence])
    2838               infeasibilities[iSequence]=COIN_INDEXED_REALLY_TINY_ELEMENT;
    2839           }
    2840         } else {
    2841           if (fabs(djValue) > freeTolerance) {
    2842             // we are going to bias towards free (but only if reasonable)
    2843             djValue *= FREE_BIAS;
    2844             if (!infeasibilities[iSequence])
    2845               newInfeas[numberNewInfeas++]=iSequence;
    2846             // store square in list
    2847             infeasibilities[iSequence] = djValue * djValue;
    2848           } else {
    2849             if (infeasibilities[iSequence])
    2850               infeasibilities[iSequence]=COIN_INDEXED_REALLY_TINY_ELEMENT;
    2851           }
    2852         }
    2853       }
    2854       if (fabs(tableauValue)>zeroTolerance&&updateWeights) {
    2855         double modification=0.0;
    2856         for (CoinBigIndex j = start; j < next; j++) {
    2857           int jRow = row[j];
    2858           modification += piWeights[jRow] * elementByColumn[j];
    2859         }
    2860         double thisWeight = weights[iSequence];
    2861         double pivot = tableauValue;
    2862         double pivotSquared = pivot * pivot;
    2863         thisWeight += pivotSquared * devex - pivot * modification;
    2864         if (thisWeight < DEVEX_TRY_NORM) {
    2865           if (referenceIn < 0.0) {
    2866             // steepest
    2867             thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
    2868           } else {
    2869             // exact
    2870             thisWeight = referenceIn * pivotSquared;
    2871             if (isReference(iSequence))
    2872               thisWeight += 1.0;
    2873             thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
    2874           }
    2875         }
    2876         weights[iSequence] = thisWeight;
    2877       }
    2878     }
    2879   }
    2880   spareColumn1.setTempNumElementsPartition(iBlock,numberNewInfeas);
    2881   int bestSequence=-1;
     2813        reducedCost[iSequence] = djValue;
     2814        if (iStatus < 4) {
     2815          djValue *= mult[iStatus];
     2816          if (djValue < 0.0) {
     2817            if (!infeasibilities[iSequence])
     2818              newInfeas[numberNewInfeas++] = iSequence;
     2819            infeasibilities[iSequence] = djValue * djValue;
     2820          } else {
     2821            if (infeasibilities[iSequence])
     2822              infeasibilities[iSequence] = COIN_INDEXED_REALLY_TINY_ELEMENT;
     2823          }
     2824        } else {
     2825          if (fabs(djValue) > freeTolerance) {
     2826            // we are going to bias towards free (but only if reasonable)
     2827            djValue *= FREE_BIAS;
     2828            if (!infeasibilities[iSequence])
     2829              newInfeas[numberNewInfeas++] = iSequence;
     2830            // store square in list
     2831            infeasibilities[iSequence] = djValue * djValue;
     2832          } else {
     2833            if (infeasibilities[iSequence])
     2834              infeasibilities[iSequence] = COIN_INDEXED_REALLY_TINY_ELEMENT;
     2835          }
     2836        }
     2837      }
     2838      if (fabs(tableauValue) > zeroTolerance && updateWeights) {
     2839        double modification = 0.0;
     2840        for (CoinBigIndex j = start; j < next; j++) {
     2841          int jRow = row[j];
     2842          modification += piWeights[jRow] * elementByColumn[j];
     2843        }
     2844        double thisWeight = weights[iSequence];
     2845        double pivot = tableauValue;
     2846        double pivotSquared = pivot * pivot;
     2847        thisWeight += pivotSquared * devex - pivot * modification;
     2848        if (thisWeight < DEVEX_TRY_NORM) {
     2849          if (referenceIn < 0.0) {
     2850            // steepest
     2851            thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
     2852          } else {
     2853            // exact
     2854            thisWeight = referenceIn * pivotSquared;
     2855            if (isReference(iSequence))
     2856              thisWeight += 1.0;
     2857            thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
     2858          }
     2859        }
     2860        weights[iSequence] = thisWeight;
     2861      }
     2862    }
     2863  }
     2864  spareColumn1.setTempNumElementsPartition(iBlock, numberNewInfeas);
     2865  int bestSequence = -1;
    28822866#if 0
    28832867  if (!iBlock)
     
    29002884   otherwise use updateForDjs
    29012885*/
    2902 int
    2903 AbcMatrix::primalColumnSparseDouble(int iBlock,CoinPartitionedVector & updateForTableauRow,
    2904                               CoinPartitionedVector & updateForDjs,
    2905                               const CoinIndexedVector & updateForWeights,
    2906                               CoinPartitionedVector & spareColumn1,
    2907                               double * infeasibilities,
    2908                               double referenceIn, double devex,
    2909                               // Array for exact devex to say what is in reference framework
    2910                               unsigned int * reference,
    2911                               double * weights, double scaleFactor) const
     2886int AbcMatrix::primalColumnSparseDouble(int iBlock, CoinPartitionedVector &updateForTableauRow,
     2887  CoinPartitionedVector &updateForDjs,
     2888  const CoinIndexedVector &updateForWeights,
     2889  CoinPartitionedVector &spareColumn1,
     2890  double *infeasibilities,
     2891  double referenceIn, double devex,
     2892  // Array for exact devex to say what is in reference framework
     2893  unsigned int *reference,
     2894  double *weights, double scaleFactor) const
    29122895{
    29132896  int maximumRows = model_->maximumAbcNumberRows();
    2914   int first=blockStart_[iBlock];
     2897  int first = blockStart_[iBlock];
    29152898  //int last=blockStart_[iBlock+1];
    2916   const double *  COIN_RESTRICT piTableau=updateForTableauRow.denseVector();
     2899  const double *COIN_RESTRICT piTableau = updateForTableauRow.denseVector();
    29172900  //const double *  COIN_RESTRICT pi=updateForDjs.denseVector();
    2918   const double *  COIN_RESTRICT piWeights=updateForWeights.denseVector();
    2919   const unsigned char *  COIN_RESTRICT internalStatus = model_->internalStatus();
    2920   double * COIN_RESTRICT reducedCost = model_->djRegion();
     2901  const double *COIN_RESTRICT piWeights = updateForWeights.denseVector();
     2902  const unsigned char *COIN_RESTRICT internalStatus = model_->internalStatus();
     2903  double *COIN_RESTRICT reducedCost = model_->djRegion();
    29212904  double tolerance = model_->currentDualTolerance();
    29222905  // use spareColumn to track new infeasibilities
    2923   int * COIN_RESTRICT newInfeas = spareColumn1.getIndices()+first;
    2924   int numberNewInfeas=0;
     2906  int *COIN_RESTRICT newInfeas = spareColumn1.getIndices() + first;
     2907  int numberNewInfeas = 0;
    29252908  // we can't really trust infeasibilities if there is dual error
    29262909  // this coding has to mimic coding in checkDualSolution
    29272910  double error = CoinMin(1.0e-2, model_->largestDualError());
    29282911  // allow tolerance at least slightly bigger than standard
    2929   tolerance = tolerance +  error;
    2930   double mult[2]={1.0,-1.0};
    2931   bool updateWeights=devex!=0.0;
    2932   int numberSlacks=updateForTableauRow.getNumElements();
    2933   const int *  COIN_RESTRICT piIndexTableau = updateForTableauRow.getIndices();
    2934 #define isReference(i) (((reference[i>>5] >> (i & 31)) & 1) != 0)
     2912  tolerance = tolerance + error;
     2913  double mult[2] = { 1.0, -1.0 };
     2914  bool updateWeights = devex != 0.0;
     2915  int numberSlacks = updateForTableauRow.getNumElements();
     2916  const int *COIN_RESTRICT piIndexTableau = updateForTableauRow.getIndices();
     2917#define isReference(i) (((reference[i >> 5] >> (i & 31)) & 1) != 0)
    29352918  // Scale factor is other way round so tableau row is 1.0* while dj update is scaleFactor*
    2936   assert (scaleFactor);
     2919  assert(scaleFactor);
    29372920  if (!iBlock) {
    2938     if (scaleFactor!=COIN_DBL_MAX) {
    2939       for (int i=0;i<numberSlacks;i++) {
    2940         int iRow=piIndexTableau[i];
    2941         int iStatus=internalStatus[iRow]&7;
    2942         if (iStatus<4) {
    2943           double value=reducedCost[iRow];
    2944           value += scaleFactor*piTableau[iRow]; // sign?
    2945           reducedCost[iRow]=value;
    2946           value *= mult[iStatus];
    2947           if (value<0.0) {
    2948             if (!infeasibilities[iRow])
    2949               newInfeas[numberNewInfeas++]=iRow;
    2950 #ifdef CLP_PRIMAL_SLACK_MULTIPLIER 
    2951             infeasibilities[iRow]=value*value*CLP_PRIMAL_SLACK_MULTIPLIER;
     2921    if (scaleFactor != COIN_DBL_MAX) {
     2922      for (int i = 0; i < numberSlacks; i++) {
     2923        int iRow = piIndexTableau[i];
     2924        int iStatus = internalStatus[iRow] & 7;
     2925        if (iStatus < 4) {
     2926          double value = reducedCost[iRow];
     2927          value += scaleFactor * piTableau[iRow]; // sign?
     2928          reducedCost[iRow] = value;
     2929          value *= mult[iStatus];
     2930          if (value < 0.0) {
     2931            if (!infeasibilities[iRow])
     2932              newInfeas[numberNewInfeas++] = iRow;
     2933#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
     2934            infeasibilities[iRow] = value * value * CLP_PRIMAL_SLACK_MULTIPLIER;
    29522935#else
    2953             infeasibilities[iRow]=value*value;
    2954 #endif
    2955           } else {
    2956             if (infeasibilities[iRow])
    2957               infeasibilities[iRow]=COIN_INDEXED_REALLY_TINY_ELEMENT;
    2958           }
    2959         }
     2936            infeasibilities[iRow] = value * value;
     2937#endif
     2938          } else {
     2939            if (infeasibilities[iRow])
     2940              infeasibilities[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
     2941          }
     2942        }
    29602943      }
    29612944    }
    29622945    if (updateWeights) {
    2963       for (int i=0;i<numberSlacks;i++) {
    2964         int iRow=piIndexTableau[i];
    2965         double modification=piWeights[iRow];
    2966         double thisWeight = weights[iRow];
    2967         double pivot = piTableau[iRow];
    2968         double pivotSquared = pivot * pivot;
    2969         thisWeight += pivotSquared * devex - pivot * modification;
    2970         if (thisWeight < DEVEX_TRY_NORM) {
    2971           if (referenceIn < 0.0) {
    2972             // steepest
    2973             thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
    2974           } else {
    2975             // exact
    2976             thisWeight = referenceIn * pivotSquared;
    2977             if (isReference(iRow))
    2978               thisWeight += 1.0;
    2979             thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
    2980           }
    2981         }
    2982         weights[iRow] = thisWeight;
    2983       }
    2984     }
    2985     first=maximumRows;
     2946      for (int i = 0; i < numberSlacks; i++) {
     2947        int iRow = piIndexTableau[i];
     2948        double modification = piWeights[iRow];
     2949        double thisWeight = weights[iRow];
     2950        double pivot = piTableau[iRow];
     2951        double pivotSquared = pivot * pivot;
     2952        thisWeight += pivotSquared * devex - pivot * modification;
     2953        if (thisWeight < DEVEX_TRY_NORM) {
     2954          if (referenceIn < 0.0) {
     2955            // steepest
     2956            thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
     2957          } else {
     2958            // exact
     2959            thisWeight = referenceIn * pivotSquared;
     2960            if (isReference(iRow))
     2961              thisWeight += 1.0;
     2962            thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
     2963          }
     2964        }
     2965        weights[iRow] = thisWeight;
     2966      }
     2967    }
     2968    first = maximumRows;
    29862969  }
    29872970  double zeroTolerance = model_->zeroTolerance();
    2988   double freeTolerance = FREE_ACCEPT * tolerance;;
     2971  double freeTolerance = FREE_ACCEPT * tolerance;
     2972  ;
    29892973  // get matrix data pointers
    2990   const int * COIN_RESTRICT row = matrix_->getIndices();
    2991   const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts()-maximumRows;
    2992   const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
     2974  const int *COIN_RESTRICT row = matrix_->getIndices();
     2975  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts() - maximumRows;
     2976  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
    29932977  // get tableau row
    2994   int *  COIN_RESTRICT index = updateForTableauRow.getIndices()+first;
    2995   double *  COIN_RESTRICT array = updateForTableauRow.denseVector();
    2996   int numberRows=model_->numberRows();
    2997   const CoinBigIndex * COIN_RESTRICT rowStart = rowStart_+iBlock*numberRows;
    2998   const CoinBigIndex * COIN_RESTRICT rowEnd = rowStart+numberRows;
    2999   const double * COIN_RESTRICT element = element_;
    3000   const int * COIN_RESTRICT column = column_;
    3001   int numberInRow=0;
    3002   if (numberSlacks>2) {
    3003     for (int i=0;i<numberSlacks;i++) {
    3004       int iRow=piIndexTableau[i];
     2978  int *COIN_RESTRICT index = updateForTableauRow.getIndices() + first;
     2979  double *COIN_RESTRICT array = updateForTableauRow.denseVector();
     2980  int numberRows = model_->numberRows();
     2981  const CoinBigIndex *COIN_RESTRICT rowStart = rowStart_ + iBlock * numberRows;
     2982  const CoinBigIndex *COIN_RESTRICT rowEnd = rowStart + numberRows;
     2983  const double *COIN_RESTRICT element = element_;
     2984  const int *COIN_RESTRICT column = column_;
     2985  int numberInRow = 0;
     2986  if (numberSlacks > 2) {
     2987    for (int i = 0; i < numberSlacks; i++) {
     2988      int iRow = piIndexTableau[i];
    30052989      double piValue = piTableau[iRow];
    30062990      CoinBigIndex end = rowEnd[iRow];
    3007       for (CoinBigIndex j=rowStart[iRow];j<end;j++) {
    3008         int iSequence = column[j];
    3009         double value = element[j]*piValue;
    3010         double oldValue=array[iSequence];
    3011         value += oldValue;
    3012         if (!oldValue) {
    3013           array[iSequence]= value;
    3014           index[numberInRow++]=iSequence;
    3015         } else if (value) {     
    3016           array[iSequence] = value;
    3017         } else {
    3018           array[iSequence] = COIN_INDEXED_REALLY_TINY_ELEMENT;
    3019         }
    3020       }
    3021     }
    3022   } else if (numberSlacks==2) {
     2991      for (CoinBigIndex j = rowStart[iRow]; j < end; j++) {
     2992        int iSequence = column[j];
     2993        double value = element[j] * piValue;
     2994        double oldValue = array[iSequence];
     2995        value += oldValue;
     2996        if (!oldValue) {
     2997          array[iSequence] = value;
     2998          index[numberInRow++] = iSequence;
     2999        } else if (value) {
     3000          array[iSequence] = value;
     3001        } else {
     3002          array[iSequence] = COIN_INDEXED_REALLY_TINY_ELEMENT;
     3003        }
     3004      }
     3005    }
     3006  } else if (numberSlacks == 2) {
    30233007    int iRow0 = piIndexTableau[0];
    30243008    int iRow1 = piIndexTableau[1];
    30253009    CoinBigIndex end0 = rowEnd[iRow0];
    30263010    CoinBigIndex end1 = rowEnd[iRow1];
    3027     if (end0-rowStart[iRow0]>end1-rowStart[iRow1]) {
    3028       int temp=iRow0;
    3029       iRow0=iRow1;
    3030       iRow1=temp;
     3011    if (end0 - rowStart[iRow0] > end1 - rowStart[iRow1]) {
     3012      int temp = iRow0;
     3013      iRow0 = iRow1;
     3014      iRow1 = temp;
    30313015    }
    30323016    CoinBigIndex start = rowStart[iRow0];
    30333017    CoinBigIndex end = rowEnd[iRow0];
    30343018    double piValue = piTableau[iRow0];
    3035     for (CoinBigIndex j=start;j<end;j++) {
     3019    for (CoinBigIndex j = start; j < end; j++) {
    30363020      int iSequence = column[j];
    30373021      double value = element[j];
    3038       array[iSequence]= piValue*value;
    3039       index[numberInRow++]=iSequence;
     3022      array[iSequence] = piValue * value;
     3023      index[numberInRow++] = iSequence;
    30403024    }
    30413025    start = rowStart[iRow1];
    30423026    end = rowEnd[iRow1];
    30433027    piValue = piTableau[iRow1];
    3044     for (CoinBigIndex j=start;j<end;j++) {
     3028    for (CoinBigIndex j = start; j < end; j++) {
    30453029      int iSequence = column[j];
    30463030      double value = element[j];
    30473031      value *= piValue;
    30483032      if (!array[iSequence]) {
    3049         array[iSequence]= value;
    3050         index[numberInRow++]=iSequence;
    3051       } else { 
    3052         value += array[iSequence];
    3053         if (value)
    3054           array[iSequence]= value;
    3055         else
    3056           array[iSequence] = COIN_INDEXED_REALLY_TINY_ELEMENT;
     3033        array[iSequence] = value;
     3034        index[numberInRow++] = iSequence;
     3035      } else {
     3036        value += array[iSequence];
     3037        if (value)
     3038          array[iSequence] = value;
     3039        else
     3040          array[iSequence] = COIN_INDEXED_REALLY_TINY_ELEMENT;
    30573041      }
    30583042    }
     
    30623046    CoinBigIndex end = rowEnd[iRow0];
    30633047    double piValue = piTableau[iRow0];
    3064     for (CoinBigIndex j=start;j<end;j++) {
     3048    for (CoinBigIndex j = start; j < end; j++) {
    30653049      int iSequence = column[j];
    30663050      double value = element[j];
    3067       array[iSequence]= piValue*value;
    3068       index[numberInRow++]=iSequence;
    3069     }
    3070   }
    3071   bool updateDjs=scaleFactor!=COIN_DBL_MAX;
    3072   for (int iLook =0 ;iLook<numberInRow;iLook++) {
    3073     int iSequence=index[iLook];
    3074     unsigned char iStatus=internalStatus[iSequence]&7;
    3075     if (iStatus<6) {
     3051      array[iSequence] = piValue * value;
     3052      index[numberInRow++] = iSequence;
     3053    }
     3054  }
     3055  bool updateDjs = scaleFactor != COIN_DBL_MAX;
     3056  for (int iLook = 0; iLook < numberInRow; iLook++) {
     3057    int iSequence = index[iLook];
     3058    unsigned char iStatus = internalStatus[iSequence] & 7;
     3059    if (iStatus < 6) {
    30763060      double tableauValue = array[iSequence];
    3077       array[iSequence]=0.0;
    3078       double djValue=reducedCost[iSequence];
    3079       djValue += tableauValue*scaleFactor; // sign?
     3061      array[iSequence] = 0.0;
     3062      double djValue = reducedCost[iSequence];
     3063      djValue += tableauValue * scaleFactor; // sign?
    30803064      if (updateDjs) {
    3081         reducedCost[iSequence]=djValue;
    3082         if (iStatus<4) {
    3083           djValue *= mult[iStatus];
    3084           if (djValue<0.0) {
    3085             if (!infeasibilities[iSequence])
    3086               newInfeas[numberNewInfeas++]=iSequence;
    3087             infeasibilities[iSequence]=djValue*djValue;
    3088           } else {
    3089             if (infeasibilities[iSequence])
    3090               infeasibilities[iSequence]=COIN_INDEXED_REALLY_TINY_ELEMENT;
    3091           }
    3092         } else {
    3093           if (fabs(djValue) > freeTolerance) {
    3094             // we are going to bias towards free (but only if reasonable)
    3095             djValue *= FREE_BIAS;
    3096             if (!infeasibilities[iSequence])
    3097               newInfeas[numberNewInfeas++]=iSequence;
    3098             // store square in list
    3099             infeasibilities[iSequence] = djValue * djValue;
    3100           } else {
    3101             if (infeasibilities[iSequence])
    3102               infeasibilities[iSequence]=COIN_INDEXED_REALLY_TINY_ELEMENT;
    3103           }
    3104         }
    3105       }
    3106       if (fabs(tableauValue)>zeroTolerance&&updateWeights) {
    3107         CoinBigIndex start = columnStart[iSequence];
    3108         CoinBigIndex next = columnStart[iSequence+1];
    3109         double modification=0.0;
    3110         for (CoinBigIndex j = start; j < next; j++) {
    3111           int jRow = row[j];
    3112           modification += piWeights[jRow] * elementByColumn[j];
    3113         }
    3114         double thisWeight = weights[iSequence];
    3115         double pivot = tableauValue;
    3116         double pivotSquared = pivot * pivot;
    3117         thisWeight += pivotSquared * devex - pivot * modification;
    3118         if (thisWeight < DEVEX_TRY_NORM) {
    3119           if (referenceIn < 0.0) {
    3120             // steepest
    3121             thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
    3122           } else {
    3123             // exact
    3124             thisWeight = referenceIn * pivotSquared;
    3125             if (isReference(iSequence))
    3126               thisWeight += 1.0;
    3127             thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
    3128           }
    3129         }
    3130         weights[iSequence] = thisWeight;
     3065        reducedCost[iSequence] = djValue;
     3066        if (iStatus < 4) {
     3067          djValue *= mult[iStatus];
     3068          if (djValue < 0.0) {
     3069            if (!infeasibilities[iSequence])
     3070              newInfeas[numberNewInfeas++] = iSequence;
     3071            infeasibilities[iSequence] = djValue * djValue;
     3072          } else {
     3073            if (infeasibilities[iSequence])
     3074              infeasibilities[iSequence] = COIN_INDEXED_REALLY_TINY_ELEMENT;
     3075          }
     3076        } else {
     3077          if (fabs(djValue) > freeTolerance) {
     3078            // we are going to bias towards free (but only if reasonable)
     3079            djValue *= FREE_BIAS;
     3080            if (!infeasibilities[iSequence])
     3081              newInfeas[numberNewInfeas++] = iSequence;
     3082            // store square in list
     3083            infeasibilities[iSequence] = djValue * djValue;
     3084          } else {
     3085            if (infeasibilities[iSequence])
     3086              infeasibilities[iSequence] = COIN_INDEXED_REALLY_TINY_ELEMENT;
     3087          }
     3088        }
     3089      }
     3090      if (fabs(tableauValue) > zeroTolerance && updateWeights) {
     3091        CoinBigIndex start = columnStart[iSequence];
     3092        CoinBigIndex next = columnStart[iSequence + 1];
     3093        double modification = 0.0;
     3094        for (CoinBigIndex j = start; j < next; j++) {
     3095          int jRow = row[j];
     3096          modification += piWeights[jRow] * elementByColumn[j];
     3097        }
     3098        double thisWeight = weights[iSequence];
     3099        double pivot = tableauValue;
     3100        double pivotSquared = pivot * pivot;
     3101        thisWeight += pivotSquared * devex - pivot * modification;
     3102        if (thisWeight < DEVEX_TRY_NORM) {
     3103          if (referenceIn < 0.0) {
     3104            // steepest
     3105            thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
     3106          } else {
     3107            // exact
     3108            thisWeight = referenceIn * pivotSquared;
     3109            if (isReference(iSequence))
     3110              thisWeight += 1.0;
     3111            thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
     3112          }
     3113        }
     3114        weights[iSequence] = thisWeight;
    31313115      }
    31323116    } else {
    3133       array[iSequence]=0.0;
    3134     }
    3135   }
    3136   spareColumn1.setTempNumElementsPartition(iBlock,numberNewInfeas);
    3137   int bestSequence=-1;
     3117      array[iSequence] = 0.0;
     3118    }
     3119  }
     3120  spareColumn1.setTempNumElementsPartition(iBlock, numberNewInfeas);
     3121  int bestSequence = -1;
    31383122#if 0
    31393123  if (!iBlock)
     
    31673151   Returns best
    31683152*/
    3169 int
    3170 AbcMatrix::primalColumnDouble(CoinPartitionedVector & updateForTableauRow,
    3171                               CoinPartitionedVector & updateForDjs,
    3172                               const CoinIndexedVector & updateForWeights,
    3173                               CoinPartitionedVector & spareColumn1,
    3174                               CoinIndexedVector & infeasible,
    3175                               double referenceIn, double devex,
    3176                               // Array for exact devex to say what is in reference framework
    3177                               unsigned int * reference,
    3178                               double * weights, double scaleFactor) const
     3153int AbcMatrix::primalColumnDouble(CoinPartitionedVector &updateForTableauRow,
     3154  CoinPartitionedVector &updateForDjs,
     3155  const CoinIndexedVector &updateForWeights,
     3156  CoinPartitionedVector &spareColumn1,
     3157  CoinIndexedVector &infeasible,
     3158  double referenceIn, double devex,
     3159  // Array for exact devex to say what is in reference framework
     3160  unsigned int *reference,
     3161  double *weights, double scaleFactor) const
    31793162{
    31803163  //int maximumRows = model_->maximumAbcNumberRows();
    3181    // rebalance
    3182    rebalance();
     3164  // rebalance
     3165  rebalance();
    31833166#ifdef PRICE_IN_ABC_MATRIX
    3184    int which[NUMBER_BLOCKS];
    3185 #endif
    3186    double * infeasibilities=infeasible.denseVector();
    3187    int bestSequence=-1;
    3188    // see if worth using row copy
    3189    int maximumRows=model_->maximumAbcNumberRows();
    3190    int number=updateForTableauRow.getNumElements();
     3167  int which[NUMBER_BLOCKS];
     3168#endif
     3169  double *infeasibilities = infeasible.denseVector();
     3170  int bestSequence = -1;
     3171  // see if worth using row copy
     3172  int maximumRows = model_->maximumAbcNumberRows();
     3173  int number = updateForTableauRow.getNumElements();
    31913174#ifdef GCC_4_9
    3192    assert (number);
     3175  assert(number);
    31933176#else
    3194    if (!number) {
    3195      printf("Null tableau row!\n");
    3196    }
    3197 #endif
    3198    bool useRowCopy = (gotRowCopy()&&(number<<2)<maximumRows);
    3199    //useRowCopy=false;
    3200    if (!scaleFactor)
    3201      useRowCopy=false; // look again later
    3202    const int * starts;
    3203    int numberBlocks;
    3204    if (useRowCopy) {
    3205      starts=blockStart_;
    3206      numberBlocks=numberRowBlocks_;
    3207    } else {
    3208      starts=startColumnBlock_;
    3209      numberBlocks=numberColumnBlocks_;
    3210    }
    3211    if (useRowCopy) {
    3212      for (int i=0;i<numberBlocks;i++)
     3177  if (!number) {
     3178    printf("Null tableau row!\n");
     3179  }
     3180#endif
     3181  bool useRowCopy = (gotRowCopy() && (number << 2) < maximumRows);
     3182  //useRowCopy=false;
     3183  if (!scaleFactor)
     3184    useRowCopy = false; // look again later
     3185  const int *starts;
     3186  int numberBlocks;
     3187  if (useRowCopy) {
     3188    starts = blockStart_;
     3189    numberBlocks = numberRowBlocks_;
     3190  } else {
     3191    starts = startColumnBlock_;
     3192    numberBlocks = numberColumnBlocks_;
     3193  }
     3194  if (useRowCopy) {
     3195    for (int i = 0; i < numberBlocks; i++)
    32133196#ifdef PRICE_IN_ABC_MATRIX
    3214        which[i]=
    3215 #endif
    3216          cilk_spawn primalColumnSparseDouble(i,updateForTableauRow,updateForDjs,updateForWeights,
    3217                                        spareColumn1,
    3218                                        infeasibilities,referenceIn,devex,reference,weights,scaleFactor);
    3219      cilk_sync;
    3220    } else {
    3221      for (int i=0;i<numberBlocks;i++)
     3197      which[i] =
     3198#endif
     3199        cilk_spawn primalColumnSparseDouble(i, updateForTableauRow, updateForDjs, updateForWeights,
     3200          spareColumn1,
     3201          infeasibilities, referenceIn, devex, reference, weights, scaleFactor);
     3202    cilk_sync;
     3203  } else {
     3204    for (int i = 0; i < numberBlocks; i++)
    32223205#ifdef PRICE_IN_ABC_MATRIX
    3223        which[i]=
    3224 #endif
    3225          cilk_spawn primalColumnDouble(i,updateForTableauRow,updateForDjs,updateForWeights,
    3226                                        spareColumn1,
    3227                                        infeasibilities,referenceIn,devex,reference,weights,scaleFactor);
    3228      cilk_sync;
    3229    }
     3206      which[i] =
     3207#endif
     3208        cilk_spawn primalColumnDouble(i, updateForTableauRow, updateForDjs, updateForWeights,
     3209          spareColumn1,
     3210          infeasibilities, referenceIn, devex, reference, weights, scaleFactor);
     3211    cilk_sync;
     3212  }
    32303213#ifdef PRICE_IN_ABC_MATRIX
    3231    double bestValue=model_->dualTolerance();
    3232    int sequenceIn[8]={-1,-1,-1,-1,-1,-1,-1,-1};
    3233 #endif
    3234    assert (numberColumnBlocks_<=8);
     3214  double bestValue = model_->dualTolerance();
     3215  int sequenceIn[8] = { -1, -1, -1, -1, -1, -1, -1, -1 };
     3216#endif
     3217  assert(numberColumnBlocks_ <= 8);
    32353218#ifdef PRICE_IN_ABC_MATRIX
    3236    double weightedDj[8];
    3237    int nPossible=0;
    3238 #endif
    3239    //const double * reducedCost = model_->djRegion();
    3240    // use spareColumn to track new infeasibilities
    3241    int numberInfeas=infeasible.getNumElements();
    3242    int * COIN_RESTRICT infeas = infeasible.getIndices();
    3243    const int * COIN_RESTRICT newInfeasAll = spareColumn1.getIndices();
    3244    for (int i=0;i<numberColumnBlocks_;i++) {
    3245      const int * COIN_RESTRICT newInfeas = newInfeasAll+starts[i];
    3246      int numberNewInfeas=spareColumn1.getNumElements(i);
    3247      spareColumn1.setTempNumElementsPartition(i,0);
    3248      for (int j=0;j<numberNewInfeas;j++)
    3249        infeas[numberInfeas++]=newInfeas[j];
     3219  double weightedDj[8];
     3220  int nPossible = 0;
     3221#endif
     3222  //const double * reducedCost = model_->djRegion();
     3223  // use spareColumn to track new infeasibilities
     3224  int numberInfeas = infeasible.getNumElements();
     3225  int *COIN_RESTRICT infeas = infeasible.getIndices();
     3226  const int *COIN_RESTRICT newInfeasAll = spareColumn1.getIndices();
     3227  for (int i = 0; i < numberColumnBlocks_; i++) {
     3228    const int *COIN_RESTRICT newInfeas = newInfeasAll + starts[i];
     3229    int numberNewInfeas = spareColumn1.getNumElements(i);
     3230    spareColumn1.setTempNumElementsPartition(i, 0);
     3231    for (int j = 0; j < numberNewInfeas; j++)
     3232      infeas[numberInfeas++] = newInfeas[j];
    32503233#ifdef PRICE_IN_ABC_MATRIX
    3251      if (which[i]>=0) {
    3252        int iSequence=which[i];
    3253        double value=fabs(infeasibilities[iSequence]/weights[iSequence]);
    3254        if (value>bestValue) {
    3255          bestValue=value;
    3256          bestSequence=which[i];
    3257        }
    3258        weightedDj[nPossible]=-value;
    3259        sequenceIn[nPossible++]=iSequence;
    3260      }
    3261 #endif
    3262    }
    3263    infeasible.setNumElements(numberInfeas);
     3234    if (which[i] >= 0) {
     3235      int iSequence = which[i];
     3236      double value = fabs(infeasibilities[iSequence] / weights[iSequence]);
     3237      if (value > bestValue) {
     3238        bestValue = value;
     3239        bestSequence = which[i];
     3240      }
     3241      weightedDj[nPossible] = -value;
     3242      sequenceIn[nPossible++] = iSequence;
     3243    }
     3244#endif
     3245  }
     3246  infeasible.setNumElements(numberInfeas);
    32643247#ifdef PRICE_IN_ABC_MATRIX
    3265    CoinSort_2(weightedDj,weightedDj+nPossible,sequenceIn);
    3266    //model_->setMultipleSequenceIn(sequenceIn);
    3267 #endif
    3268    return bestSequence;
    3269 
     3248  CoinSort_2(weightedDj, weightedDj + nPossible, sequenceIn);
     3249  //model_->setMultipleSequenceIn(sequenceIn);
     3250#endif
     3251  return bestSequence;
    32703252}
    32713253// Partial pricing
    3272 void
    3273 AbcMatrix::partialPricing(double startFraction, double endFraction,
    3274                           int & bestSequence, int & numberWanted)
     3254void AbcMatrix::partialPricing(double startFraction, double endFraction,
     3255  int &bestSequence, int &numberWanted)
    32753256{
    32763257  numberWanted = currentWanted_;
    32773258  int maximumRows = model_->maximumAbcNumberRows();
    32783259  // get matrix data pointers
    3279   const int * COIN_RESTRICT row = matrix_->getIndices();
    3280   const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts()-maximumRows;
    3281   const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
    3282   int numberColumns=model_->numberColumns();
    3283   int start = static_cast<int> (startFraction * numberColumns);
    3284   int end = CoinMin(static_cast<int> (endFraction * numberColumns + 1), numberColumns);
     3260  const int *COIN_RESTRICT row = matrix_->getIndices();
     3261  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts() - maximumRows;
     3262  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
     3263  int numberColumns = model_->numberColumns();
     3264  int start = static_cast< int >(startFraction * numberColumns);
     3265  int end = CoinMin(static_cast< int >(endFraction * numberColumns + 1), numberColumns);
    32853266  // adjust
    32863267  start += maximumRows;
    32873268  end += maximumRows;
    32883269  double tolerance = model_->currentDualTolerance();
    3289   double * reducedCost = model_->djRegion();
    3290   const double * duals = model_->dualRowSolution();
    3291   const double * cost = model_->costRegion();
     3270  double *reducedCost = model_->djRegion();
     3271  const double *duals = model_->dualRowSolution();
     3272  const double *cost = model_->costRegion();
    32923273  double bestDj;
    32933274  if (bestSequence >= 0)
     
    33033284      double value;
    33043285      AbcSimplex::Status status = model_->getInternalStatus(iSequence);
    3305       switch(status) {
     3286      switch (status) {
    33063287      case AbcSimplex::basic:
    33073288      case AbcSimplex::isFixed:
    3308         break;
     3289        break;
    33093290      case AbcSimplex::isFree:
    33103291      case AbcSimplex::superBasic:
    3311         value = cost[iSequence];
    3312         for (CoinBigIndex j = columnStart[iSequence];
    3313              j < columnStart[iSequence+1]; j++) {
    3314           int jRow = row[j];
    3315           value -= duals[jRow] * elementByColumn[j];
    3316         }
    3317         value = fabs(value);
    3318         if (value > FREE_ACCEPT * tolerance) {
    3319           numberWanted--;
    3320           // we are going to bias towards free (but only if reasonable)
    3321           value *= FREE_BIAS;
    3322           if (value > bestDj) {
    3323             // check flagged variable and correct dj
    3324             if (!model_->flagged(iSequence)) {
    3325               bestDj = value;
    3326               bestSequence = iSequence;
    3327             } else {
    3328               // just to make sure we don't exit before got something
    3329               numberWanted++;
    3330             }
    3331           }
    3332         }
    3333         break;
     3292        value = cost[iSequence];
     3293        for (CoinBigIndex j = columnStart[iSequence];
     3294             j < columnStart[iSequence + 1]; j++) {
     3295          int jRow = row[j];
     3296          value -= duals[jRow] * elementByColumn[j];
     3297        }
     3298        value = fabs(value);
     3299        if (value > FREE_ACCEPT * tolerance) {
     3300          numberWanted--;
     3301          // we are going to bias towards free (but only if reasonable)
     3302          value *= FREE_BIAS;
     3303          if (value > bestDj) {
     3304            // check flagged variable and correct dj
     3305            if (!model_->flagged(iSequence)) {
     3306              bestDj = value;
     3307              bestSequence = iSequence;
     3308            } else {
     3309              // just to make sure we don't exit before got something
     3310              numberWanted++;
     3311            }
     3312          }
     3313        }
     3314        break;
    33343315      case AbcSimplex::atUpperBound:
    3335         value = cost[iSequence];
    3336         // scaled
    3337         for (CoinBigIndex j = columnStart[iSequence];
    3338              j < columnStart[iSequence+1]; j++) {
    3339           int jRow = row[j];
    3340           value -= duals[jRow] * elementByColumn[j];
    3341         }
    3342         if (value > tolerance) {
    3343           numberWanted--;
    3344           if (value > bestDj) {
    3345             // check flagged variable and correct dj
    3346             if (!model_->flagged(iSequence)) {
    3347               bestDj = value;
    3348               bestSequence = iSequence;
    3349             } else {
    3350               // just to make sure we don't exit before got something
    3351               numberWanted++;
    3352             }
    3353           }
    3354         }
    3355         break;
     3316        value = cost[iSequence];
     3317        // scaled
     3318        for (CoinBigIndex j = columnStart[iSequence];
     3319             j < columnStart[iSequence + 1]; j++) {
     3320          int jRow = row[j];
     3321          value -= duals[jRow] * elementByColumn[j];
     3322        }
     3323        if (value > tolerance) {
     3324          numberWanted--;
     3325          if (value > bestDj) {
     3326            // check flagged variable and correct dj
     3327            if (!model_->flagged(iSequence)) {
     3328              bestDj = value;
     3329              bestSequence = iSequence;
     3330            } else {
     3331              // just to make sure we don't exit before got something
     3332              numberWanted++;
     3333            }
     3334          }
     3335        }
     3336        break;
    33563337      case AbcSimplex::atLowerBound:
    3357         value = cost[iSequence];
    3358         for (CoinBigIndex j = columnStart[iSequence];
    3359              j < columnStart[iSequence+1]; j++) {
    3360           int jRow = row[j];
    3361           value -= duals[jRow] * elementByColumn[j];
    3362         }
    3363         value = -value;
    3364         if (value > tolerance) {
    3365           numberWanted--;
    3366           if (value > bestDj) {
    3367             // check flagged variable and correct dj
    3368             if (!model_->flagged(iSequence)) {
    3369               bestDj = value;
    3370               bestSequence = iSequence;
    3371             } else {
    3372               // just to make sure we don't exit before got something
    3373               numberWanted++;
    3374             }
    3375           }
    3376         }
    3377         break;
     3338        value = cost[iSequence];
     3339        for (CoinBigIndex j = columnStart[iSequence];
     3340             j < columnStart[iSequence + 1]; j++) {
     3341          int jRow = row[j];
     3342          value -= duals[jRow] * elementByColumn[j];
     3343        }
     3344        value = -value;
     3345        if (value > tolerance) {
     3346          numberWanted--;
     3347          if (value > bestDj) {
     3348            // check flagged variable and correct dj
     3349            if (!model_->flagged(iSequence)) {
     3350              bestDj = value;
     3351              bestSequence = iSequence;
     3352            } else {
     3353              // just to make sure we don't exit before got something
     3354              numberWanted++;
     3355            }
     3356          }
     3357        }
     3358        break;
    33783359      }
    33793360    }
     
    33893370    double value = cost[bestSequence];
    33903371    for (CoinBigIndex j = columnStart[bestSequence];
    3391          j < columnStart[bestSequence+1]; j++) {
     3372         j < columnStart[bestSequence + 1]; j++) {
    33923373      int jRow = row[j];
    33933374      value -= duals[jRow] * elementByColumn[j];
     
    34023383   just for indices Already in z.
    34033384   Note - z always packed mode */
    3404 void
    3405 AbcMatrix::subsetTransposeTimes(const CoinIndexedVector & x,
    3406                                 CoinIndexedVector & z) const
     3385void AbcMatrix::subsetTransposeTimes(const CoinIndexedVector &x,
     3386  CoinIndexedVector &z) const
    34073387{
    34083388  int maximumRows = model_->maximumAbcNumberRows();
    34093389  // get matrix data pointers
    3410   const int * COIN_RESTRICT row = matrix_->getIndices();
    3411   const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts()-maximumRows;
    3412   const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
    3413   const double *  COIN_RESTRICT other=x.denseVector();
     3390  const int *COIN_RESTRICT row = matrix_->getIndices();
     3391  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts() - maximumRows;
     3392  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
     3393  const double *COIN_RESTRICT other = x.denseVector();
    34143394  int numberNonZero = z.getNumElements();
    3415   int *  COIN_RESTRICT index = z.getIndices();
    3416   double *  COIN_RESTRICT array = z.denseVector();
    3417   int numberRows=model_->numberRows();
    3418   for (int i=0;i<numberNonZero;i++) {
    3419     int iSequence=index[i];
    3420     if (iSequence>=numberRows) {
     3395  int *COIN_RESTRICT index = z.getIndices();
     3396  double *COIN_RESTRICT array = z.denseVector();
     3397  int numberRows = model_->numberRows();
     3398  for (int i = 0; i < numberNonZero; i++) {
     3399    int iSequence = index[i];
     3400    if (iSequence >= numberRows) {
    34213401      CoinBigIndex start = columnStart[iSequence];
    3422       CoinBigIndex next = columnStart[iSequence+1];
     3402      CoinBigIndex next = columnStart[iSequence + 1];
    34233403      double value = 0.0;
    34243404      for (CoinBigIndex j = start; j < next; j++) {
    3425         int jRow = row[j];
    3426         value -= other[jRow] * elementByColumn[j];
    3427       }
    3428       array[i]=value;
    3429     } else {
    3430       array[i]=-other[iSequence];
    3431     }
    3432   }
    3433 }
    3434 // Return <code>-x *A</code> in <code>z</code>
    3435 void
    3436 AbcMatrix::transposeTimes(const CoinIndexedVector & x,
    3437                                 CoinIndexedVector & z) const
     3405        int jRow = row[j];
     3406        value -= other[jRow] * elementByColumn[j];
     3407      }
     3408      array[i] = value;
     3409    } else {
     3410      array[i] = -other[iSequence];
     3411    }
     3412  }
     3413}
     3414// Return <code>-x *A</code> in <code>z</code>
     3415void AbcMatrix::transposeTimes(const CoinIndexedVector &x,
     3416  CoinIndexedVector &z) const
    34383417{
    34393418  int maximumRows = model_->maximumAbcNumberRows();
    34403419  // get matrix data pointers
    3441   const int * COIN_RESTRICT row = matrix_->getIndices();
    3442   const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts()-maximumRows;
    3443   const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
    3444   const double *  COIN_RESTRICT other=x.denseVector();
     3420  const int *COIN_RESTRICT row = matrix_->getIndices();
     3421  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts() - maximumRows;
     3422  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
     3423  const double *COIN_RESTRICT other = x.denseVector();
    34453424  int numberNonZero = 0;
    3446   int *  COIN_RESTRICT index = z.getIndices();
    3447   double *  COIN_RESTRICT array = z.denseVector();
    3448   int numberColumns=model_->numberColumns();
    3449   double zeroTolerance=model_->zeroTolerance();
    3450   for (int iSequence=maximumRows;iSequence<maximumRows+numberColumns;iSequence++) {
     3425  int *COIN_RESTRICT index = z.getIndices();
     3426  double *COIN_RESTRICT array = z.denseVector();
     3427  int numberColumns = model_->numberColumns();
     3428  double zeroTolerance = model_->zeroTolerance();
     3429  for (int iSequence = maximumRows; iSequence < maximumRows + numberColumns; iSequence++) {
    34513430    CoinBigIndex start = columnStart[iSequence];
    3452     CoinBigIndex next = columnStart[iSequence+1];
     3431    CoinBigIndex next = columnStart[iSequence + 1];
    34533432    double value = 0.0;
    34543433    for (CoinBigIndex j = start; j < next; j++) {
     
    34563435      value -= other[jRow] * elementByColumn[j];
    34573436    }
    3458     if (fabs(value)>zeroTolerance) {
     3437    if (fabs(value) > zeroTolerance) {
    34593438      // TEMP
    3460       array[iSequence-maximumRows]=value;
    3461       index[numberNonZero++]=iSequence-maximumRows;
     3439