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

    r2367 r2385  
    33// Corporation and others.  All Rights Reserved.
    44// This code is licensed under the terms of the Eclipse Public License (EPL).
    5 
    6 
    75
    86#include <cstdio>
     
    3331#ifdef COIN_PREFETCH
    3432#if 1
    35 #define coin_prefetch(mem) \
    36          __asm__ __volatile__ ("prefetchnta %0" : : "m" (*(reinterpret_cast<char *>(mem))))
    37 #define coin_prefetch_const(mem) \
    38          __asm__ __volatile__ ("prefetchnta %0" : : "m" (*(reinterpret_cast<const char *>(mem))))
     33#define coin_prefetch(mem)              \
     34  __asm__ __volatile__("prefetchnta %0" \
     35                       :                \
     36                       : "m"(*(reinterpret_cast< char * >(mem))))
     37#define coin_prefetch_const(mem)        \
     38  __asm__ __volatile__("prefetchnta %0" \
     39                       :                \
     40                       : "m"(*(reinterpret_cast< const char * >(mem))))
    3941#else
    40 #define coin_prefetch(mem) \
    41          __asm__ __volatile__ ("prefetch %0" : : "m" (*(reinterpret_cast<char *>(mem))))
    42 #define coin_prefetch_const(mem) \
    43          __asm__ __volatile__ ("prefetch %0" : : "m" (*(reinterpret_cast<const char *>(mem))))
     42#define coin_prefetch(mem)           \
     43  __asm__ __volatile__("prefetch %0" \
     44                       :             \
     45                       : "m"(*(reinterpret_cast< char * >(mem))))
     46#define coin_prefetch_const(mem)     \
     47  __asm__ __volatile__("prefetch %0" \
     48                       :             \
     49                       : "m"(*(reinterpret_cast< const char * >(mem))))
    4450#endif
    4551#else
     
    5258// Constructors / Destructor / Assignment
    5359//#############################################################################
    54 static void debug3(int iColumn,double &thisWeight,double pivotSquared,double devex,
    55                    double pivot,double modification,double oldWeight)
     60static void debug3(int iColumn, double &thisWeight, double pivotSquared, double devex,
     61  double pivot, double modification, double oldWeight)
    5662{
    5763  //double newWeight=oldWeight+pivotSquared * devex + pivot * modification;
     
    5965  //printf("YY %d %.30g %.30g %.30g %.30g %.30g %.30g\n",
    6066  //       iColumn,thisWeight,pivotSquared,devex,pivot,modification,oldWeight);
    61   thisWeight=oldWeight;
     67  thisWeight = oldWeight;
    6268  thisWeight += pivotSquared * devex + pivot * modification;
    6369}
     
    6571// Default Constructor
    6672//-------------------------------------------------------------------
    67 ClpPackedMatrix::ClpPackedMatrix ()
    68      : ClpMatrixBase(),
    69        matrix_(NULL),
    70        numberActiveColumns_(0),
    71        flags_(2),
    72        rowCopy_(NULL),
    73       columnCopy_(NULL)
    74 {
    75      setType(1);
     73ClpPackedMatrix::ClpPackedMatrix()
     74  : ClpMatrixBase()
     75  , matrix_(NULL)
     76  , numberActiveColumns_(0)
     77  , flags_(2)
     78  , rowCopy_(NULL)
     79  , columnCopy_(NULL)
     80{
     81  setType(1);
    7682}
    7783
     
    7985// Copy constructor
    8086//-------------------------------------------------------------------
    81 ClpPackedMatrix::ClpPackedMatrix (const ClpPackedMatrix & rhs)
    82      : ClpMatrixBase(rhs)
     87ClpPackedMatrix::ClpPackedMatrix(const ClpPackedMatrix &rhs)
     88  : ClpMatrixBase(rhs)
    8389{
    8490#ifdef DO_CHECK_FLAGS
    85      rhs.checkFlags(0);
     91  rhs.checkFlags(0);
    8692#endif
    8793#ifndef COIN_SPARSE_MATRIX
    88      // Guaranteed no gaps or small elements
    89      matrix_ = new CoinPackedMatrix(*(rhs.matrix_),-1,0) ;
    90      flags_ = rhs.flags_&(~0x02) ;
     94  // Guaranteed no gaps or small elements
     95  matrix_ = new CoinPackedMatrix(*(rhs.matrix_), -1, 0);
     96  flags_ = rhs.flags_ & (~0x02);
    9197#else
    92      // Gaps & small elements preserved
    93      matrix_ = new CoinPackedMatrix(*(rhs.matrix_),0,0) ;
    94      flags_ = rhs.flags_ ;
    95      if (matrix_->hasGaps()) flags_ |= 0x02 ;
    96 #endif
    97      numberActiveColumns_ = rhs.numberActiveColumns_;
    98      int numberRows = matrix_->getNumRows();
    99      if (rhs.rhsOffset_ && numberRows) {
    100           rhsOffset_ = ClpCopyOfArray(rhs.rhsOffset_, numberRows);
    101      } else {
    102           rhsOffset_ = NULL;
    103      }
    104      if (rhs.rowCopy_) {
    105           assert ((flags_ & 4) != 0);
    106           rowCopy_ = new ClpPackedMatrix2(*rhs.rowCopy_);
    107      } else {
    108           rowCopy_ = NULL;
    109      }
    110      if (rhs.columnCopy_) {
    111           assert ((flags_&(8 + 16)) == 8 + 16);
    112           columnCopy_ = new ClpPackedMatrix3(*rhs.columnCopy_);
    113      } else {
    114           columnCopy_ = NULL;
    115      }
     98  // Gaps & small elements preserved
     99  matrix_ = new CoinPackedMatrix(*(rhs.matrix_), 0, 0);
     100  flags_ = rhs.flags_;
     101  if (matrix_->hasGaps())
     102    flags_ |= 0x02;
     103#endif
     104  numberActiveColumns_ = rhs.numberActiveColumns_;
     105  int numberRows = matrix_->getNumRows();
     106  if (rhs.rhsOffset_ && numberRows) {
     107    rhsOffset_ = ClpCopyOfArray(rhs.rhsOffset_, numberRows);
     108  } else {
     109    rhsOffset_ = NULL;
     110  }
     111  if (rhs.rowCopy_) {
     112    assert((flags_ & 4) != 0);
     113    rowCopy_ = new ClpPackedMatrix2(*rhs.rowCopy_);
     114  } else {
     115    rowCopy_ = NULL;
     116  }
     117  if (rhs.columnCopy_) {
     118    assert((flags_ & (8 + 16)) == 8 + 16);
     119    columnCopy_ = new ClpPackedMatrix3(*rhs.columnCopy_);
     120  } else {
     121    columnCopy_ = NULL;
     122  }
    116123#ifdef DO_CHECK_FLAGS
    117      checkFlags(0);
     124  checkFlags(0);
    118125#endif
    119126}
     
    121128// assign matrix (for space reasons)
    122129//-------------------------------------------------------------------
    123 ClpPackedMatrix::ClpPackedMatrix (CoinPackedMatrix * rhs)
    124      : ClpMatrixBase()
    125 {
    126      matrix_ = rhs;
    127      flags_ = ((matrix_->hasGaps())?0x02:0) ;
    128      numberActiveColumns_ = matrix_->getNumCols();
    129      rowCopy_ = NULL;
    130      columnCopy_ = NULL;
    131      setType(1);
     130ClpPackedMatrix::ClpPackedMatrix(CoinPackedMatrix *rhs)
     131  : ClpMatrixBase()
     132{
     133  matrix_ = rhs;
     134  flags_ = ((matrix_->hasGaps()) ? 0x02 : 0);
     135  numberActiveColumns_ = matrix_->getNumCols();
     136  rowCopy_ = NULL;
     137  columnCopy_ = NULL;
     138  setType(1);
    132139#ifdef DO_CHECK_FLAGS
    133      checkFlags(0);
    134 #endif
     140  checkFlags(0);
     141#endif
     142}
    135143
    136 }
    137 
    138 ClpPackedMatrix::ClpPackedMatrix (const CoinPackedMatrix & rhs)
    139      : ClpMatrixBase()
     144ClpPackedMatrix::ClpPackedMatrix(const CoinPackedMatrix &rhs)
     145  : ClpMatrixBase()
    140146{
    141147#ifndef COIN_SPARSE_MATRIX
    142      matrix_ = new CoinPackedMatrix(rhs,-1,0) ;
    143      flags_ = 0 ;
     148  matrix_ = new CoinPackedMatrix(rhs, -1, 0);
     149  flags_ = 0;
    144150#else
    145      matrix_ = new CoinPackedMatrix(rhs,0,0) ;
    146      flags_ = ((matrix_->hasGaps())?0x02:0) ;
    147 #endif
    148      numberActiveColumns_ = matrix_->getNumCols();
    149      rowCopy_ = NULL;
    150      columnCopy_ = NULL;
    151      setType(1);
     151  matrix_ = new CoinPackedMatrix(rhs, 0, 0);
     152  flags_ = ((matrix_->hasGaps()) ? 0x02 : 0);
     153#endif
     154  numberActiveColumns_ = matrix_->getNumCols();
     155  rowCopy_ = NULL;
     156  columnCopy_ = NULL;
     157  setType(1);
    152158#ifdef DO_CHECK_FLAGS
    153      checkFlags(0);
    154 #endif
    155 
     159  checkFlags(0);
     160#endif
    156161}
    157162
     
    159164// Destructor
    160165//-------------------------------------------------------------------
    161 ClpPackedMatrix::~ClpPackedMatrix ()
    162 {
    163      delete matrix_;
    164      delete rowCopy_;
    165      delete columnCopy_;
     166ClpPackedMatrix::~ClpPackedMatrix()
     167{
     168  delete matrix_;
     169  delete rowCopy_;
     170  delete columnCopy_;
    166171}
    167172
     
    170175//-------------------------------------------------------------------
    171176ClpPackedMatrix &
    172 ClpPackedMatrix::operator=(const ClpPackedMatrix& rhs)
    173 {
    174      if (this != &rhs) {
    175           ClpMatrixBase::operator=(rhs);
    176           delete matrix_;
     177ClpPackedMatrix::operator=(const ClpPackedMatrix &rhs)
     178{
     179  if (this != &rhs) {
     180    ClpMatrixBase::operator=(rhs);
     181    delete matrix_;
    177182#ifndef COIN_SPARSE_MATRIX
    178           matrix_ = new CoinPackedMatrix(*(rhs.matrix_),-1,0) ;
    179           flags_ = rhs.flags_&(~0x02) ;
     183    matrix_ = new CoinPackedMatrix(*(rhs.matrix_), -1, 0);
     184    flags_ = rhs.flags_ & (~0x02);
    180185#else
    181           matrix_ = new CoinPackedMatrix(*(rhs.matrix_),0,0) ;
    182           flags_ = rhs.flags_ ;
    183           if (matrix_->hasGaps()) flags_ |= 0x02 ;
    184 #endif
    185           numberActiveColumns_ = rhs.numberActiveColumns_;
    186           delete rowCopy_;
    187           delete columnCopy_;
    188           if (rhs.rowCopy_) {
    189                assert ((flags_ & 4) != 0);
    190                rowCopy_ = new ClpPackedMatrix2(*rhs.rowCopy_);
    191           } else {
    192                rowCopy_ = NULL;
    193           }
    194           if (rhs.columnCopy_) {
    195                assert ((flags_&(8 + 16)) == 8 + 16);
    196                columnCopy_ = new ClpPackedMatrix3(*rhs.columnCopy_);
    197           } else {
    198                columnCopy_ = NULL;
    199           }
     186    matrix_ = new CoinPackedMatrix(*(rhs.matrix_), 0, 0);
     187    flags_ = rhs.flags_;
     188    if (matrix_->hasGaps())
     189      flags_ |= 0x02;
     190#endif
     191    numberActiveColumns_ = rhs.numberActiveColumns_;
     192    delete rowCopy_;
     193    delete columnCopy_;
     194    if (rhs.rowCopy_) {
     195      assert((flags_ & 4) != 0);
     196      rowCopy_ = new ClpPackedMatrix2(*rhs.rowCopy_);
     197    } else {
     198      rowCopy_ = NULL;
     199    }
     200    if (rhs.columnCopy_) {
     201      assert((flags_ & (8 + 16)) == 8 + 16);
     202      columnCopy_ = new ClpPackedMatrix3(*rhs.columnCopy_);
     203    } else {
     204      columnCopy_ = NULL;
     205    }
    200206#ifdef DO_CHECK_FLAGS
    201           checkFlags(0);
    202 #endif
    203      }
    204      return *this;
     207    checkFlags(0);
     208#endif
     209  }
     210  return *this;
    205211}
    206212//-------------------------------------------------------------------
    207213// Clone
    208214//-------------------------------------------------------------------
    209 ClpMatrixBase * ClpPackedMatrix::clone() const
    210 {
    211      return new ClpPackedMatrix(*this);
     215ClpMatrixBase *ClpPackedMatrix::clone() const
     216{
     217  return new ClpPackedMatrix(*this);
    212218}
    213219// Copy contents - resizing if necessary - otherwise re-use memory
    214 void
    215 ClpPackedMatrix::copy(const ClpPackedMatrix * rhs)
    216 {
    217      //*this = *rhs;
    218      assert(numberActiveColumns_ == rhs->numberActiveColumns_);
    219      assert (matrix_->isColOrdered() == rhs->matrix_->isColOrdered());
    220      matrix_->copyReuseArrays(*rhs->matrix_);
     220void ClpPackedMatrix::copy(const ClpPackedMatrix *rhs)
     221{
     222  //*this = *rhs;
     223  assert(numberActiveColumns_ == rhs->numberActiveColumns_);
     224  assert(matrix_->isColOrdered() == rhs->matrix_->isColOrdered());
     225  matrix_->copyReuseArrays(*rhs->matrix_);
    221226#ifdef DO_CHECK_FLAGS
    222      checkFlags(0);
     227  checkFlags(0);
    223228#endif
    224229}
     
    226231   and order is as given */
    227232ClpMatrixBase *
    228 ClpPackedMatrix::subsetClone (int numberRows, const int * whichRows,
    229                               int numberColumns,
    230                               const int * whichColumns) const
    231 {
    232      return new ClpPackedMatrix(*this, numberRows, whichRows,
    233                                 numberColumns, whichColumns);
     233ClpPackedMatrix::subsetClone(int numberRows, const int *whichRows,
     234  int numberColumns,
     235  const int *whichColumns) const
     236{
     237  return new ClpPackedMatrix(*this, numberRows, whichRows,
     238    numberColumns, whichColumns);
    234239}
    235240/* Subset constructor (without gaps).  Duplicates are allowed
    236241   and order is as given */
    237 ClpPackedMatrix::ClpPackedMatrix (
    238      const ClpPackedMatrix & rhs,
    239      int numberRows, const int * whichRows,
    240      int numberColumns, const int * whichColumns)
    241      : ClpMatrixBase(rhs)
    242 {
    243      matrix_ = new CoinPackedMatrix(*(rhs.matrix_), numberRows, whichRows,
    244                                     numberColumns, whichColumns);
    245      numberActiveColumns_ = matrix_->getNumCols();
    246      rowCopy_ = NULL;
    247      flags_ = rhs.flags_&(~0x02) ; // no gaps
    248      columnCopy_ = NULL;
     242ClpPackedMatrix::ClpPackedMatrix(
     243  const ClpPackedMatrix &rhs,
     244  int numberRows, const int *whichRows,
     245  int numberColumns, const int *whichColumns)
     246  : ClpMatrixBase(rhs)
     247{
     248  matrix_ = new CoinPackedMatrix(*(rhs.matrix_), numberRows, whichRows,
     249    numberColumns, whichColumns);
     250  numberActiveColumns_ = matrix_->getNumCols();
     251  rowCopy_ = NULL;
     252  flags_ = rhs.flags_ & (~0x02); // no gaps
     253  columnCopy_ = NULL;
    249254#ifdef DO_CHECK_FLAGS
    250      checkFlags(0);
    251 #endif
    252 }
    253 ClpPackedMatrix::ClpPackedMatrix (
    254      const CoinPackedMatrix & rhs,
    255      int numberRows, const int * whichRows,
    256      int numberColumns, const int * whichColumns)
    257      : ClpMatrixBase()
    258 {
    259      matrix_ = new CoinPackedMatrix(rhs, numberRows, whichRows,
    260                                     numberColumns, whichColumns);
    261      numberActiveColumns_ = matrix_->getNumCols();
    262      rowCopy_ = NULL;
    263      flags_ = 0 ; // no gaps
    264      columnCopy_ = NULL;
    265      setType(1);
     255  checkFlags(0);
     256#endif
     257}
     258ClpPackedMatrix::ClpPackedMatrix(
     259  const CoinPackedMatrix &rhs,
     260  int numberRows, const int *whichRows,
     261  int numberColumns, const int *whichColumns)
     262  : ClpMatrixBase()
     263{
     264  matrix_ = new CoinPackedMatrix(rhs, numberRows, whichRows,
     265    numberColumns, whichColumns);
     266  numberActiveColumns_ = matrix_->getNumCols();
     267  rowCopy_ = NULL;
     268  flags_ = 0; // no gaps
     269  columnCopy_ = NULL;
     270  setType(1);
    266271#ifdef DO_CHECK_FLAGS
    267      checkFlags(0);
     272  checkFlags(0);
    268273#endif
    269274}
     
    273278ClpPackedMatrix::reverseOrderedCopy() const
    274279{
    275      ClpPackedMatrix * copy = new ClpPackedMatrix();
    276      copy->matrix_ = new CoinPackedMatrix();
    277      copy->matrix_->setExtraGap(0.0);
    278      copy->matrix_->setExtraMajor(0.0);
    279      copy->matrix_->reverseOrderedCopyOf(*matrix_);
    280      //copy->matrix_->removeGaps();
    281      copy->numberActiveColumns_ = copy->matrix_->getNumCols();
    282      copy->flags_ = flags_&(~0x02) ; // no gaps
     280  ClpPackedMatrix *copy = new ClpPackedMatrix();
     281  copy->matrix_ = new CoinPackedMatrix();
     282  copy->matrix_->setExtraGap(0.0);
     283  copy->matrix_->setExtraMajor(0.0);
     284  copy->matrix_->reverseOrderedCopyOf(*matrix_);
     285  //copy->matrix_->removeGaps();
     286  copy->numberActiveColumns_ = copy->matrix_->getNumCols();
     287  copy->flags_ = flags_ & (~0x02); // no gaps
    283288#ifdef DO_CHECK_FLAGS
    284      checkFlags(0);
    285 #endif
    286      return copy;
     289  checkFlags(0);
     290#endif
     291  return copy;
    287292}
    288293//unscaled versions
    289 void
    290 ClpPackedMatrix::times(double scalar,
    291                        const double * COIN_RESTRICT x, double * COIN_RESTRICT y) const
    292 {
    293      int iRow, iColumn;
    294      // get matrix data pointers
    295      const int * COIN_RESTRICT row = matrix_->getIndices();
    296      const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
    297      const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
    298      //memset(y,0,matrix_->getNumRows()*sizeof(double));
    299      assert (((flags_&0x02) != 0) == matrix_->hasGaps());
    300      if (!(flags_ & 2)) {
    301           for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
    302                CoinBigIndex j;
    303                double value = x[iColumn];
    304                if (value) {
    305                     CoinBigIndex start = columnStart[iColumn];
    306                     CoinBigIndex end = columnStart[iColumn+1];
    307                     value *= scalar;
    308                     for (j = start; j < end; j++) {
    309                          iRow = row[j];
    310                          y[iRow] += value * elementByColumn[j];
    311                     }
    312                }
    313           }
    314      } else {
    315           const int * columnLength = matrix_->getVectorLengths();
    316           for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
    317                CoinBigIndex j;
    318                double value = x[iColumn];
    319                if (value) {
    320                     CoinBigIndex start = columnStart[iColumn];
    321                     CoinBigIndex end = start + columnLength[iColumn];
    322                     value *= scalar;
    323                     for (j = start; j < end; j++) {
    324                          iRow = row[j];
    325                          y[iRow] += value * elementByColumn[j];
    326                     }
    327                }
    328           }
    329      }
     294void ClpPackedMatrix::times(double scalar,
     295  const double *COIN_RESTRICT x, double *COIN_RESTRICT y) const
     296{
     297  int iRow, iColumn;
     298  // get matrix data pointers
     299  const int *COIN_RESTRICT row = matrix_->getIndices();
     300  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     301  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
     302  //memset(y,0,matrix_->getNumRows()*sizeof(double));
     303  assert(((flags_ & 0x02) != 0) == matrix_->hasGaps());
     304  if (!(flags_ & 2)) {
     305    for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
     306      CoinBigIndex j;
     307      double value = x[iColumn];
     308      if (value) {
     309        CoinBigIndex start = columnStart[iColumn];
     310        CoinBigIndex end = columnStart[iColumn + 1];
     311        value *= scalar;
     312        for (j = start; j < end; j++) {
     313          iRow = row[j];
     314          y[iRow] += value * elementByColumn[j];
     315        }
     316      }
     317    }
     318  } else {
     319    const int *columnLength = matrix_->getVectorLengths();
     320    for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
     321      CoinBigIndex j;
     322      double value = x[iColumn];
     323      if (value) {
     324        CoinBigIndex start = columnStart[iColumn];
     325        CoinBigIndex end = start + columnLength[iColumn];
     326        value *= scalar;
     327        for (j = start; j < end; j++) {
     328          iRow = row[j];
     329          y[iRow] += value * elementByColumn[j];
     330        }
     331      }
     332    }
     333  }
    330334}
    331335#if ABOCA_LITE
    332336static void
    333 transposeTimesBit(clpTempInfo & info)
    334 {
    335   const CoinBigIndex * COIN_RESTRICT columnStart=info.start;
    336   const int * COIN_RESTRICT row = info.row;
    337   const double * COIN_RESTRICT elementByColumn = info.element;
    338   double * COIN_RESTRICT y = info.spare;
    339   const double * COIN_RESTRICT x = info.work;
     337transposeTimesBit(clpTempInfo &info)
     338{
     339  const CoinBigIndex *COIN_RESTRICT columnStart = info.start;
     340  const int *COIN_RESTRICT row = info.row;
     341  const double *COIN_RESTRICT elementByColumn = info.element;
     342  double *COIN_RESTRICT y = info.spare;
     343  const double *COIN_RESTRICT x = info.work;
    340344  int first = info.startColumn;
    341   int last = first+info.numberToDo;
     345  int last = first + info.numberToDo;
    342346  CoinBigIndex start = columnStart[first];
    343347  for (int iColumn = first; iColumn < last; iColumn++) {
    344348    CoinBigIndex j;
    345     CoinBigIndex next = columnStart[iColumn+1];
     349    CoinBigIndex next = columnStart[iColumn + 1];
    346350    double value = y[iColumn];
    347351    for (j = start; j < next; j++) {
     
    354358}
    355359#endif
    356 void
    357 ClpPackedMatrix::transposeTimes(double scalar,
    358                                 const double * COIN_RESTRICT x, double * COIN_RESTRICT y) const
    359 {
    360      int iColumn;
    361      // get matrix data pointers
    362      const int * COIN_RESTRICT row = matrix_->getIndices();
    363      const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
    364      const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
    365      if (!(flags_ & 2)) {
    366           if (scalar == -1.0) {
     360void ClpPackedMatrix::transposeTimes(double scalar,
     361  const double *COIN_RESTRICT x, double *COIN_RESTRICT y) const
     362{
     363  int iColumn;
     364  // get matrix data pointers
     365  const int *COIN_RESTRICT row = matrix_->getIndices();
     366  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     367  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
     368  if (!(flags_ & 2)) {
     369    if (scalar == -1.0) {
    367370#if ABOCA_LITE
    368             int numberThreads=abcState();
    369             if (numberThreads) {
    370             clpTempInfo info[ABOCA_LITE];
    371             int chunk = (numberActiveColumns_+numberThreads-1)/numberThreads;
    372             int n=0;
    373             for (int i=0;i<numberThreads;i++) {
    374               info[i].spare=y;
    375               info[i].work=const_cast<double *>(x);
    376               info[i].startColumn=n;
    377               info[i].element=elementByColumn;
    378               info[i].start=columnStart;
    379               info[i].row=row;
    380               info[i].numberToDo=CoinMin(chunk,numberActiveColumns_-n);
    381               n += chunk;
    382             }
    383             for (int i=0;i<numberThreads;i++) {
    384               cilk_spawn transposeTimesBit(info[i]);
    385             }
    386             cilk_sync;
    387             } else {
    388 #endif
    389                CoinBigIndex start = columnStart[0];
    390                for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
    391                     CoinBigIndex j;
    392                     CoinBigIndex next = columnStart[iColumn+1];
    393                     double value = y[iColumn];
    394                     for (j = start; j < next; j++) {
    395                          int jRow = row[j];
    396                          value -= x[jRow] * elementByColumn[j];
    397                     }
    398                     start = next;
    399                     y[iColumn] = value;
    400                }
     371      int numberThreads = abcState();
     372      if (numberThreads) {
     373        clpTempInfo info[ABOCA_LITE];
     374        int chunk = (numberActiveColumns_ + numberThreads - 1) / numberThreads;
     375        int n = 0;
     376        for (int i = 0; i < numberThreads; i++) {
     377          info[i].spare = y;
     378          info[i].work = const_cast< double * >(x);
     379          info[i].startColumn = n;
     380          info[i].element = elementByColumn;
     381          info[i].start = columnStart;
     382          info[i].row = row;
     383          info[i].numberToDo = CoinMin(chunk, numberActiveColumns_ - n);
     384          n += chunk;
     385        }
     386        for (int i = 0; i < numberThreads; i++) {
     387          cilk_spawn transposeTimesBit(info[i]);
     388        }
     389        cilk_sync;
     390      } else {
     391#endif
     392        CoinBigIndex start = columnStart[0];
     393        for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
     394          CoinBigIndex j;
     395          CoinBigIndex next = columnStart[iColumn + 1];
     396          double value = y[iColumn];
     397          for (j = start; j < next; j++) {
     398            int jRow = row[j];
     399            value -= x[jRow] * elementByColumn[j];
     400          }
     401          start = next;
     402          y[iColumn] = value;
     403        }
    401404#if ABOCA_LITE
    402             }
    403 #endif
    404           } else {
    405                CoinBigIndex start = columnStart[0];
    406                for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
    407                     CoinBigIndex j;
    408                     CoinBigIndex next = columnStart[iColumn+1];
    409                     double value = 0.0;
    410                     for (j = start; j < next; j++) {
    411                          int jRow = row[j];
    412                          value += x[jRow] * elementByColumn[j];
    413                     }
    414                     start = next;
    415                     y[iColumn] += value * scalar;
    416                }
    417           }
    418      } else {
    419           const int * columnLength = matrix_->getVectorLengths();
     405      }
     406#endif
     407    } else {
     408      CoinBigIndex start = columnStart[0];
     409      for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
     410        CoinBigIndex j;
     411        CoinBigIndex next = columnStart[iColumn + 1];
     412        double value = 0.0;
     413        for (j = start; j < next; j++) {
     414          int jRow = row[j];
     415          value += x[jRow] * elementByColumn[j];
     416        }
     417        start = next;
     418        y[iColumn] += value * scalar;
     419      }
     420    }
     421  } else {
     422    const int *columnLength = matrix_->getVectorLengths();
     423    for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
     424      CoinBigIndex j;
     425      double value = 0.0;
     426      CoinBigIndex start = columnStart[iColumn];
     427      CoinBigIndex end = start + columnLength[iColumn];
     428      for (j = start; j < end; j++) {
     429        int jRow = row[j];
     430        value += x[jRow] * elementByColumn[j];
     431      }
     432      y[iColumn] += value * scalar;
     433    }
     434  }
     435}
     436void ClpPackedMatrix::times(double scalar,
     437  const double *COIN_RESTRICT x, double *COIN_RESTRICT y,
     438  const double *COIN_RESTRICT rowScale,
     439  const double *COIN_RESTRICT columnScale) const
     440{
     441  if (rowScale) {
     442    int iRow, iColumn;
     443    // get matrix data pointers
     444    const int *COIN_RESTRICT row = matrix_->getIndices();
     445    const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     446    const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
     447    if (!(flags_ & 2)) {
     448      for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
     449        double value = x[iColumn];
     450        if (value) {
     451          // scaled
     452          value *= scalar * columnScale[iColumn];
     453          CoinBigIndex start = columnStart[iColumn];
     454          CoinBigIndex end = columnStart[iColumn + 1];
     455          CoinBigIndex j;
     456          for (j = start; j < end; j++) {
     457            iRow = row[j];
     458            y[iRow] += value * elementByColumn[j] * rowScale[iRow];
     459          }
     460        }
     461      }
     462    } else {
     463      const int *COIN_RESTRICT columnLength = matrix_->getVectorLengths();
     464      for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
     465        double value = x[iColumn];
     466        if (value) {
     467          // scaled
     468          value *= scalar * columnScale[iColumn];
     469          CoinBigIndex start = columnStart[iColumn];
     470          CoinBigIndex end = start + columnLength[iColumn];
     471          CoinBigIndex j;
     472          for (j = start; j < end; j++) {
     473            iRow = row[j];
     474            y[iRow] += value * elementByColumn[j] * rowScale[iRow];
     475          }
     476        }
     477      }
     478    }
     479  } else {
     480    times(scalar, x, y);
     481  }
     482}
     483void ClpPackedMatrix::transposeTimes(double scalar,
     484  const double *COIN_RESTRICT x, double *COIN_RESTRICT y,
     485  const double *COIN_RESTRICT rowScale,
     486  const double *COIN_RESTRICT columnScale,
     487  double *COIN_RESTRICT spare) const
     488{
     489  if (rowScale) {
     490    int iColumn;
     491    // get matrix data pointers
     492    const int *COIN_RESTRICT row = matrix_->getIndices();
     493    const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     494    const int *COIN_RESTRICT columnLength = matrix_->getVectorLengths();
     495    const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
     496    if (!spare) {
     497      if (!(flags_ & 2)) {
     498        CoinBigIndex start = columnStart[0];
     499        if (scalar == -1.0) {
    420500          for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
    421                CoinBigIndex j;
    422                double value = 0.0;
    423                CoinBigIndex start = columnStart[iColumn];
    424                CoinBigIndex end = start + columnLength[iColumn];
    425                for (j = start; j < end; j++) {
    426                     int jRow = row[j];
    427                     value += x[jRow] * elementByColumn[j];
    428                }
    429                y[iColumn] += value * scalar;
    430           }
    431      }
    432 }
    433 void
    434 ClpPackedMatrix::times(double scalar,
    435                        const double * COIN_RESTRICT x, double * COIN_RESTRICT y,
    436                        const double * COIN_RESTRICT rowScale,
    437                        const double * COIN_RESTRICT columnScale) const
    438 {
    439      if (rowScale) {
    440           int iRow, iColumn;
    441           // get matrix data pointers
    442           const int * COIN_RESTRICT row = matrix_->getIndices();
    443           const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
    444           const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
    445           if (!(flags_ & 2)) {
    446                for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
    447                     double value = x[iColumn];
    448                     if (value) {
    449                          // scaled
    450                          value *= scalar * columnScale[iColumn];
    451                          CoinBigIndex start = columnStart[iColumn];
    452                          CoinBigIndex end = columnStart[iColumn+1];
    453                          CoinBigIndex j;
    454                          for (j = start; j < end; j++) {
    455                               iRow = row[j];
    456                               y[iRow] += value * elementByColumn[j] * rowScale[iRow];
    457                          }
    458                     }
    459                }
    460           } else {
    461                const int * COIN_RESTRICT columnLength = matrix_->getVectorLengths();
    462                for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
    463                     double value = x[iColumn];
    464                     if (value) {
    465                          // scaled
    466                          value *= scalar * columnScale[iColumn];
    467                          CoinBigIndex start = columnStart[iColumn];
    468                          CoinBigIndex end = start + columnLength[iColumn];
    469                          CoinBigIndex j;
    470                          for (j = start; j < end; j++) {
    471                               iRow = row[j];
    472                               y[iRow] += value * elementByColumn[j] * rowScale[iRow];
    473                          }
    474                     }
    475                }
    476           }
    477      } else {
    478           times(scalar, x, y);
    479      }
    480 }
    481 void
    482 ClpPackedMatrix::transposeTimes( double scalar,
    483                                  const double * COIN_RESTRICT x, double * COIN_RESTRICT y,
    484                                  const double * COIN_RESTRICT rowScale,
    485                                  const double * COIN_RESTRICT columnScale,
    486                                  double * COIN_RESTRICT spare) const
    487 {
    488      if (rowScale) {
    489           int iColumn;
    490           // get matrix data pointers
    491           const int * COIN_RESTRICT row = matrix_->getIndices();
    492           const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
    493           const int * COIN_RESTRICT columnLength = matrix_->getVectorLengths();
    494           const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
    495           if (!spare) {
    496                if (!(flags_ & 2)) {
    497                     CoinBigIndex start = columnStart[0];
    498                     if (scalar == -1.0) {
    499                          for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
    500                               CoinBigIndex j;
    501                               CoinBigIndex next = columnStart[iColumn+1];
    502                               double value = 0.0;
    503                               // scaled
    504                               for (j = start; j < next; j++) {
    505                                    int jRow = row[j];
    506                                    value += x[jRow] * elementByColumn[j] * rowScale[jRow];
    507                               }
    508                               start = next;
    509                               y[iColumn] -= value * columnScale[iColumn];
    510                          }
    511                     } else {
    512                          for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
    513                               CoinBigIndex j;
    514                               CoinBigIndex next = columnStart[iColumn+1];
    515                               double value = 0.0;
    516                               // scaled
    517                               for (j = start; j < next; j++) {
    518                                    int jRow = row[j];
    519                                    value += x[jRow] * elementByColumn[j] * rowScale[jRow];
    520                               }
    521                               start = next;
    522                               y[iColumn] += value * scalar * columnScale[iColumn];
    523                          }
    524                     }
    525                } else {
    526                     for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
    527                          CoinBigIndex j;
    528                          double value = 0.0;
    529                          // scaled
    530                          for (j = columnStart[iColumn];
    531                                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    532                               int jRow = row[j];
    533                               value += x[jRow] * elementByColumn[j] * rowScale[jRow];
    534                          }
    535                          y[iColumn] += value * scalar * columnScale[iColumn];
    536                     }
    537                }
    538           } else {
    539                // can use spare region
    540                int iRow;
    541                int numberRows = matrix_->getNumRows();
    542                for (iRow = 0; iRow < numberRows; iRow++) {
    543                     double value = x[iRow];
    544                     if (value)
    545                          spare[iRow] = value * rowScale[iRow];
    546                     else
    547                          spare[iRow] = 0.0;
    548                }
    549                if (!(flags_ & 2)) {
    550                     CoinBigIndex start = columnStart[0];
    551                     for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
    552                          CoinBigIndex j;
    553                          CoinBigIndex next = columnStart[iColumn+1];
    554                          double value = 0.0;
    555                          // scaled
    556                          for (j = start; j < next; j++) {
    557                               int jRow = row[j];
    558                               value += spare[jRow] * elementByColumn[j];
    559                          }
    560                          start = next;
    561                          y[iColumn] += value * scalar * columnScale[iColumn];
    562                     }
    563                } else {
    564                     for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
    565                          CoinBigIndex j;
    566                          double value = 0.0;
    567                          // scaled
    568                          for (j = columnStart[iColumn];
    569                                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    570                               int jRow = row[j];
    571                               value += spare[jRow] * elementByColumn[j];
    572                          }
    573                          y[iColumn] += value * scalar * columnScale[iColumn];
    574                     }
    575                }
    576                // no need to zero out
    577                //for (iRow=0;iRow<numberRows;iRow++)
    578                //spare[iRow] = 0.0;
    579           }
    580      } else {
    581           transposeTimes(scalar, x, y);
    582      }
     501            CoinBigIndex j;
     502            CoinBigIndex next = columnStart[iColumn + 1];
     503            double value = 0.0;
     504            // scaled
     505            for (j = start; j < next; j++) {
     506              int jRow = row[j];
     507              value += x[jRow] * elementByColumn[j] * rowScale[jRow];
     508            }
     509            start = next;
     510            y[iColumn] -= value * columnScale[iColumn];
     511          }
     512        } else {
     513          for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
     514            CoinBigIndex j;
     515            CoinBigIndex next = columnStart[iColumn + 1];
     516            double value = 0.0;
     517            // scaled
     518            for (j = start; j < next; j++) {
     519              int jRow = row[j];
     520              value += x[jRow] * elementByColumn[j] * rowScale[jRow];
     521            }
     522            start = next;
     523            y[iColumn] += value * scalar * columnScale[iColumn];
     524          }
     525        }
     526      } else {
     527        for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
     528          CoinBigIndex j;
     529          double value = 0.0;
     530          // scaled
     531          for (j = columnStart[iColumn];
     532               j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     533            int jRow = row[j];
     534            value += x[jRow] * elementByColumn[j] * rowScale[jRow];
     535          }
     536          y[iColumn] += value * scalar * columnScale[iColumn];
     537        }
     538      }
     539    } else {
     540      // can use spare region
     541      int iRow;
     542      int numberRows = matrix_->getNumRows();
     543      for (iRow = 0; iRow < numberRows; iRow++) {
     544        double value = x[iRow];
     545        if (value)
     546          spare[iRow] = value * rowScale[iRow];
     547        else
     548          spare[iRow] = 0.0;
     549      }
     550      if (!(flags_ & 2)) {
     551        CoinBigIndex start = columnStart[0];
     552        for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
     553          CoinBigIndex j;
     554          CoinBigIndex next = columnStart[iColumn + 1];
     555          double value = 0.0;
     556          // scaled
     557          for (j = start; j < next; j++) {
     558            int jRow = row[j];
     559            value += spare[jRow] * elementByColumn[j];
     560          }
     561          start = next;
     562          y[iColumn] += value * scalar * columnScale[iColumn];
     563        }
     564      } else {
     565        for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
     566          CoinBigIndex j;
     567          double value = 0.0;
     568          // scaled
     569          for (j = columnStart[iColumn];
     570               j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     571            int jRow = row[j];
     572            value += spare[jRow] * elementByColumn[j];
     573          }
     574          y[iColumn] += value * scalar * columnScale[iColumn];
     575        }
     576      }
     577      // no need to zero out
     578      //for (iRow=0;iRow<numberRows;iRow++)
     579      //spare[iRow] = 0.0;
     580    }
     581  } else {
     582    transposeTimes(scalar, x, y);
     583  }
    583584}
    584585#if ABOCA_LITE
    585586static void
    586 transposeTimesSubsetBit(clpTempInfo & info)
    587 {
    588   const CoinBigIndex * COIN_RESTRICT columnStart=info.start;
    589   const int * COIN_RESTRICT row = info.row;
    590   const double * COIN_RESTRICT elementByColumn = info.element;
    591   double * COIN_RESTRICT y = info.spare;
    592   const double * COIN_RESTRICT x = info.work;
    593   const int * COIN_RESTRICT which = info.which;
     587transposeTimesSubsetBit(clpTempInfo &info)
     588{
     589  const CoinBigIndex *COIN_RESTRICT columnStart = info.start;
     590  const int *COIN_RESTRICT row = info.row;
     591  const double *COIN_RESTRICT elementByColumn = info.element;
     592  double *COIN_RESTRICT y = info.spare;
     593  const double *COIN_RESTRICT x = info.work;
     594  const int *COIN_RESTRICT which = info.which;
    594595  int first = info.startColumn;
    595   int last = first+info.numberToDo;
     596  int last = first + info.numberToDo;
    596597  for (int i = first; i < last; i++) {
    597598    int iColumn = which[i];
    598599    CoinBigIndex j;
    599600    CoinBigIndex start = columnStart[iColumn];
    600     CoinBigIndex next = columnStart[iColumn+1];
     601    CoinBigIndex next = columnStart[iColumn + 1];
    601602    double value = 0.0;
    602603    for (j = start; j < next; j++) {
     
    609610}
    610611#endif
    611 void
    612 ClpPackedMatrix::transposeTimesSubset( int number,
    613                                        const int * which,
    614                                        const double * COIN_RESTRICT x, double * COIN_RESTRICT y,
    615                                        const double * COIN_RESTRICT rowScale,
    616                                        const double * COIN_RESTRICT columnScale,
    617                                        double * COIN_RESTRICT spare) const
    618 {
    619      // get matrix data pointers
    620      const int * COIN_RESTRICT row = matrix_->getIndices();
    621      const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
    622      const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
    623      if (!spare || !rowScale) {
    624           if (rowScale) {
    625                for (int jColumn = 0; jColumn < number; jColumn++) {
    626                     int iColumn = which[jColumn];
    627                     CoinBigIndex j;
    628                     CoinBigIndex start = columnStart[iColumn];
    629                     CoinBigIndex next = columnStart[iColumn+1];
    630                     double value = 0.0;
    631                     for (j = start; j < next; j++) {
    632                          int jRow = row[j];
    633                          value += x[jRow] * elementByColumn[j] * rowScale[jRow];
    634                     }
    635                     y[iColumn] -= value * columnScale[iColumn];
    636                }
    637           } else {
     612void ClpPackedMatrix::transposeTimesSubset(int number,
     613  const int *which,
     614  const double *COIN_RESTRICT x, double *COIN_RESTRICT y,
     615  const double *COIN_RESTRICT rowScale,
     616  const double *COIN_RESTRICT columnScale,
     617  double *COIN_RESTRICT spare) const
     618{
     619  // get matrix data pointers
     620  const int *COIN_RESTRICT row = matrix_->getIndices();
     621  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     622  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
     623  if (!spare || !rowScale) {
     624    if (rowScale) {
     625      for (int jColumn = 0; jColumn < number; jColumn++) {
     626        int iColumn = which[jColumn];
     627        CoinBigIndex j;
     628        CoinBigIndex start = columnStart[iColumn];
     629        CoinBigIndex next = columnStart[iColumn + 1];
     630        double value = 0.0;
     631        for (j = start; j < next; j++) {
     632          int jRow = row[j];
     633          value += x[jRow] * elementByColumn[j] * rowScale[jRow];
     634        }
     635        y[iColumn] -= value * columnScale[iColumn];
     636      }
     637    } else {
    638638#if ABOCA_LITE
    639                int numberThreads=abcState();
    640                if (numberThreads) {
    641                clpTempInfo info[ABOCA_LITE];
    642                int chunk = (number+numberThreads-1)/numberThreads;
    643                int n=0;
    644                for (int i=0;i<numberThreads;i++) {
    645                  info[i].spare=y;
    646                  info[i].work=const_cast<double *>(x);
    647                  info[i].which=const_cast<int *>(which);
    648                  info[i].startColumn=n;
    649                  info[i].element=elementByColumn;
    650                  info[i].start=columnStart;
    651                  info[i].row=row;
    652                  info[i].numberToDo=CoinMin(chunk,number-n);
    653                 n += chunk;
    654                }
    655                for (int i=0;i<numberThreads;i++) {
    656                 cilk_spawn transposeTimesSubsetBit(info[i]);
    657                }
    658                cilk_sync;
    659                } else {
    660 #endif
    661                for (int jColumn = 0; jColumn < number; jColumn++) {
    662                     int iColumn = which[jColumn];
    663                     CoinBigIndex j;
    664                     CoinBigIndex start = columnStart[iColumn];
    665                     CoinBigIndex next = columnStart[iColumn+1];
    666                     double value = 0.0;
    667                     for (j = start; j < next; j++) {
    668                          int jRow = row[j];
    669                          value += x[jRow] * elementByColumn[j];
    670                     }
    671                     y[iColumn] -= value;
    672                }
     639      int numberThreads = abcState();
     640      if (numberThreads) {
     641        clpTempInfo info[ABOCA_LITE];
     642        int chunk = (number + numberThreads - 1) / numberThreads;
     643        int n = 0;
     644        for (int i = 0; i < numberThreads; i++) {
     645          info[i].spare = y;
     646          info[i].work = const_cast< double * >(x);
     647          info[i].which = const_cast< int * >(which);
     648          info[i].startColumn = n;
     649          info[i].element = elementByColumn;
     650          info[i].start = columnStart;
     651          info[i].row = row;
     652          info[i].numberToDo = CoinMin(chunk, number - n);
     653          n += chunk;
     654        }
     655        for (int i = 0; i < numberThreads; i++) {
     656          cilk_spawn transposeTimesSubsetBit(info[i]);
     657        }
     658        cilk_sync;
     659      } else {
     660#endif
     661        for (int jColumn = 0; jColumn < number; jColumn++) {
     662          int iColumn = which[jColumn];
     663          CoinBigIndex j;
     664          CoinBigIndex start = columnStart[iColumn];
     665          CoinBigIndex next = columnStart[iColumn + 1];
     666          double value = 0.0;
     667          for (j = start; j < next; j++) {
     668            int jRow = row[j];
     669            value += x[jRow] * elementByColumn[j];
     670          }
     671          y[iColumn] -= value;
     672        }
    673673#if ABOCA_LITE
    674             }
    675 #endif
    676           }
    677      } else {
    678           // can use spare region
    679           int iRow;
    680           int numberRows = matrix_->getNumRows();
    681           for (iRow = 0; iRow < numberRows; iRow++) {
    682                double value = x[iRow];
    683                if (value)
    684                     spare[iRow] = value * rowScale[iRow];
    685                else
    686                     spare[iRow] = 0.0;
    687           }
    688           for (int jColumn = 0; jColumn < number; jColumn++) {
    689                int iColumn = which[jColumn];
    690                CoinBigIndex j;
    691                CoinBigIndex start = columnStart[iColumn];
    692                CoinBigIndex next = columnStart[iColumn+1];
    693                double value = 0.0;
    694                for (j = start; j < next; j++) {
    695                     int jRow = row[j];
    696                     value += spare[jRow] * elementByColumn[j];
    697                }
    698                y[iColumn] -= value * columnScale[iColumn];
    699           }
    700      }
     674      }
     675#endif
     676    }
     677  } else {
     678    // can use spare region
     679    int iRow;
     680    int numberRows = matrix_->getNumRows();
     681    for (iRow = 0; iRow < numberRows; iRow++) {
     682      double value = x[iRow];
     683      if (value)
     684        spare[iRow] = value * rowScale[iRow];
     685      else
     686        spare[iRow] = 0.0;
     687    }
     688    for (int jColumn = 0; jColumn < number; jColumn++) {
     689      int iColumn = which[jColumn];
     690      CoinBigIndex j;
     691      CoinBigIndex start = columnStart[iColumn];
     692      CoinBigIndex next = columnStart[iColumn + 1];
     693      double value = 0.0;
     694      for (j = start; j < next; j++) {
     695        int jRow = row[j];
     696        value += spare[jRow] * elementByColumn[j];
     697      }
     698      y[iColumn] -= value * columnScale[iColumn];
     699    }
     700  }
    701701}
    702702/* Return <code>x * A in <code>z</code>.
    703703        Squashes small elements and knows about ClpSimplex */
    704 void
    705 ClpPackedMatrix::transposeTimes(const ClpSimplex * model, double scalar,
    706                                 const CoinIndexedVector * rowArray,
    707                                 CoinIndexedVector * y,
    708                                 CoinIndexedVector * columnArray) const
    709 {
    710      columnArray->clear();
    711      double * COIN_RESTRICT pi = rowArray->denseVector();
    712      int numberNonZero = 0;
    713      int * COIN_RESTRICT index = columnArray->getIndices();
    714      double * COIN_RESTRICT array = columnArray->denseVector();
    715      int numberInRowArray = rowArray->getNumElements();
    716      // maybe I need one in OsiSimplex
    717      double zeroTolerance = model->zeroTolerance();
     704void ClpPackedMatrix::transposeTimes(const ClpSimplex *model, double scalar,
     705  const CoinIndexedVector *rowArray,
     706  CoinIndexedVector *y,
     707  CoinIndexedVector *columnArray) const
     708{
     709  columnArray->clear();
     710  double *COIN_RESTRICT pi = rowArray->denseVector();
     711  int numberNonZero = 0;
     712  int *COIN_RESTRICT index = columnArray->getIndices();
     713  double *COIN_RESTRICT array = columnArray->denseVector();
     714  int numberInRowArray = rowArray->getNumElements();
     715  // maybe I need one in OsiSimplex
     716  double zeroTolerance = model->zeroTolerance();
    718717#if 0 //def COIN_DEVELOP
    719718     if (zeroTolerance != 1.0e-13) {
     
    721720     }
    722721#endif
    723      int numberRows = model->numberRows();
    724      ClpPackedMatrix* rowCopy =
    725           static_cast< ClpPackedMatrix*>(model->rowCopy());
    726      bool packed = rowArray->packedMode();
    727      double factor = (numberRows < 100) ? 0.25 : 0.35;
    728      factor = 0.5;
    729      // We may not want to do by row if there may be cache problems
    730      // It would be nice to find L2 cache size - for moment 512K
    731      // Be slightly optimistic
    732      if (numberActiveColumns_ * sizeof(double) > 1000000) {
    733           if (numberRows * 10 < numberActiveColumns_)
    734                factor *= 0.333333333;
    735           else if (numberRows * 4 < numberActiveColumns_)
    736                factor *= 0.5;
    737           else if (numberRows * 2 < numberActiveColumns_)
    738                factor *= 0.66666666667;
    739           //if (model->numberIterations()%50==0)
    740           //printf("%d nonzero\n",numberInRowArray);
    741      }
    742      // if not packed then bias a bit more towards by column
    743      if (!packed)
    744           factor *= 0.9;
    745      assert (!y->getNumElements());
    746      double multiplierX = 0.8;
    747      double factor2 = factor * multiplierX;
    748      if (packed && rowCopy_ && numberInRowArray > 2 && numberInRowArray > factor2 * numberRows &&
    749                numberInRowArray < 0.9 * numberRows && 0) {
    750           rowCopy_->transposeTimes(model, rowCopy->matrix_, rowArray, y, columnArray);
    751           return;
    752      }
    753      if (columnCopy_)
    754        factor *= 0.7;
    755      if (numberInRowArray > factor * numberRows || !rowCopy) {
    756           // do by column
    757           // If no gaps - can do a bit faster
    758           if (!(flags_ & 2) || columnCopy_) {
    759                transposeTimesByColumn( model,  scalar,
    760                                        rowArray, y, columnArray);
    761                return;
    762           }
    763           int iColumn;
    764           // get matrix data pointers
    765           const int * COIN_RESTRICT row = matrix_->getIndices();
    766           const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
    767           const int * COIN_RESTRICT columnLength = matrix_->getVectorLengths();
    768           const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
    769           const double * COIN_RESTRICT rowScale = model->rowScale();
     722  int numberRows = model->numberRows();
     723  ClpPackedMatrix *rowCopy = static_cast< ClpPackedMatrix * >(model->rowCopy());
     724  bool packed = rowArray->packedMode();
     725  double factor = (numberRows < 100) ? 0.25 : 0.35;
     726  factor = 0.5;
     727  // We may not want to do by row if there may be cache problems
     728  // It would be nice to find L2 cache size - for moment 512K
     729  // Be slightly optimistic
     730  if (numberActiveColumns_ * sizeof(double) > 1000000) {
     731    if (numberRows * 10 < numberActiveColumns_)
     732      factor *= 0.333333333;
     733    else if (numberRows * 4 < numberActiveColumns_)
     734      factor *= 0.5;
     735    else if (numberRows * 2 < numberActiveColumns_)
     736      factor *= 0.66666666667;
     737    //if (model->numberIterations()%50==0)
     738    //printf("%d nonzero\n",numberInRowArray);
     739  }
     740  // if not packed then bias a bit more towards by column
     741  if (!packed)
     742    factor *= 0.9;
     743  assert(!y->getNumElements());
     744  double multiplierX = 0.8;
     745  double factor2 = factor * multiplierX;
     746  if (packed && rowCopy_ && numberInRowArray > 2 && numberInRowArray > factor2 * numberRows && numberInRowArray < 0.9 * numberRows && 0) {
     747    rowCopy_->transposeTimes(model, rowCopy->matrix_, rowArray, y, columnArray);
     748    return;
     749  }
     750  if (columnCopy_)
     751    factor *= 0.7;
     752  if (numberInRowArray > factor * numberRows || !rowCopy) {
     753    // do by column
     754    // If no gaps - can do a bit faster
     755    if (!(flags_ & 2) || columnCopy_) {
     756      transposeTimesByColumn(model, scalar,
     757        rowArray, y, columnArray);
     758      return;
     759    }
     760    int iColumn;
     761    // get matrix data pointers
     762    const int *COIN_RESTRICT row = matrix_->getIndices();
     763    const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     764    const int *COIN_RESTRICT columnLength = matrix_->getVectorLengths();
     765    const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
     766    const double *COIN_RESTRICT rowScale = model->rowScale();
    770767#if 0
    771768          ClpPackedMatrix * scaledMatrix = model->clpScaledMatrix();
     
    779776          }
    780777#endif
    781           if (packed) {
    782                // need to expand pi into y
    783                assert(y->capacity() >= numberRows);
    784                double * COIN_RESTRICT piOld = pi;
    785                pi = y->denseVector();
    786                const int * COIN_RESTRICT whichRow = rowArray->getIndices();
    787                int i;
    788                if (!rowScale) {
    789                     // modify pi so can collapse to one loop
    790                     if (scalar == -1.0) {
    791                          for (i = 0; i < numberInRowArray; i++) {
    792                               int iRow = whichRow[i];
    793                               pi[iRow] = -piOld[i];
    794                          }
    795                     } else {
    796                          for (i = 0; i < numberInRowArray; i++) {
    797                               int iRow = whichRow[i];
    798                               pi[iRow] = scalar * piOld[i];
    799                          }
    800                     }
    801                     for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
    802                          double value = 0.0;
    803                          CoinBigIndex j;
    804                          for (j = columnStart[iColumn];
    805                                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    806                               int iRow = row[j];
    807                               value += pi[iRow] * elementByColumn[j];
    808                          }
    809                          if (fabs(value) > zeroTolerance) {
    810                               array[numberNonZero] = value;
    811                               index[numberNonZero++] = iColumn;
    812                          }
    813                     }
    814                } else {
     778    if (packed) {
     779      // need to expand pi into y
     780      assert(y->capacity() >= numberRows);
     781      double *COIN_RESTRICT piOld = pi;
     782      pi = y->denseVector();
     783      const int *COIN_RESTRICT whichRow = rowArray->getIndices();
     784      int i;
     785      if (!rowScale) {
     786        // modify pi so can collapse to one loop
     787        if (scalar == -1.0) {
     788          for (i = 0; i < numberInRowArray; i++) {
     789            int iRow = whichRow[i];
     790            pi[iRow] = -piOld[i];
     791          }
     792        } else {
     793          for (i = 0; i < numberInRowArray; i++) {
     794            int iRow = whichRow[i];
     795            pi[iRow] = scalar * piOld[i];
     796          }
     797        }
     798        for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
     799          double value = 0.0;
     800          CoinBigIndex j;
     801          for (j = columnStart[iColumn];
     802               j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     803            int iRow = row[j];
     804            value += pi[iRow] * elementByColumn[j];
     805          }
     806          if (fabs(value) > zeroTolerance) {
     807            array[numberNonZero] = value;
     808            index[numberNonZero++] = iColumn;
     809          }
     810        }
     811      } else {
    815812#ifdef CLP_INVESTIGATE
    816                     if (model->clpScaledMatrix())
    817                          printf("scaledMatrix_ at %d of ClpPackedMatrix\n", __LINE__);
    818 #endif
    819                     // scaled
    820                     // modify pi so can collapse to one loop
    821                     if (scalar == -1.0) {
    822                          for (i = 0; i < numberInRowArray; i++) {
    823                               int iRow = whichRow[i];
    824                               pi[iRow] = -piOld[i] * rowScale[iRow];
    825                          }
    826                     } else {
    827                          for (i = 0; i < numberInRowArray; i++) {
    828                               int iRow = whichRow[i];
    829                               pi[iRow] = scalar * piOld[i] * rowScale[iRow];
    830                          }
    831                     }
    832                     for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
    833                          double value = 0.0;
    834                          CoinBigIndex j;
    835                          const double * columnScale = model->columnScale();
    836                          for (j = columnStart[iColumn];
    837                                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    838                               int iRow = row[j];
    839                               value += pi[iRow] * elementByColumn[j];
    840                          }
    841                          value *= columnScale[iColumn];
    842                          if (fabs(value) > zeroTolerance) {
    843                               array[numberNonZero] = value;
    844                               index[numberNonZero++] = iColumn;
    845                          }
    846                     }
    847                }
    848                // zero out
    849                for (i = 0; i < numberInRowArray; i++) {
    850                     int iRow = whichRow[i];
    851                     pi[iRow] = 0.0;
    852                }
    853           } else {
    854                if (!rowScale) {
    855                     if (scalar == -1.0) {
    856                          for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
    857                               double value = 0.0;
    858                               CoinBigIndex j;
    859                               for (j = columnStart[iColumn];
    860                                         j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    861                                    int iRow = row[j];
    862                                    value += pi[iRow] * elementByColumn[j];
    863                               }
    864                               if (fabs(value) > zeroTolerance) {
    865                                    index[numberNonZero++] = iColumn;
    866                                    array[iColumn] = -value;
    867                               }
    868                          }
    869                     } else {
    870                          for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
    871                               double value = 0.0;
    872                               CoinBigIndex j;
    873                               for (j = columnStart[iColumn];
    874                                         j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    875                                    int iRow = row[j];
    876                                    value += pi[iRow] * elementByColumn[j];
    877                               }
    878                               value *= scalar;
    879                               if (fabs(value) > zeroTolerance) {
    880                                    index[numberNonZero++] = iColumn;
    881                                    array[iColumn] = value;
    882                               }
    883                          }
    884                     }
    885                } else {
     813        if (model->clpScaledMatrix())
     814          printf("scaledMatrix_ at %d of ClpPackedMatrix\n", __LINE__);
     815#endif
     816        // scaled
     817        // modify pi so can collapse to one loop
     818        if (scalar == -1.0) {
     819          for (i = 0; i < numberInRowArray; i++) {
     820            int iRow = whichRow[i];
     821            pi[iRow] = -piOld[i] * rowScale[iRow];
     822          }
     823        } else {
     824          for (i = 0; i < numberInRowArray; i++) {
     825            int iRow = whichRow[i];
     826            pi[iRow] = scalar * piOld[i] * rowScale[iRow];
     827          }
     828        }
     829        for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
     830          double value = 0.0;
     831          CoinBigIndex j;
     832          const double *columnScale = model->columnScale();
     833          for (j = columnStart[iColumn];
     834               j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     835            int iRow = row[j];
     836            value += pi[iRow] * elementByColumn[j];
     837          }
     838          value *= columnScale[iColumn];
     839          if (fabs(value) > zeroTolerance) {
     840            array[numberNonZero] = value;
     841            index[numberNonZero++] = iColumn;
     842          }
     843        }
     844      }
     845      // zero out
     846      for (i = 0; i < numberInRowArray; i++) {
     847        int iRow = whichRow[i];
     848        pi[iRow] = 0.0;
     849      }
     850    } else {
     851      if (!rowScale) {
     852        if (scalar == -1.0) {
     853          for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
     854            double value = 0.0;
     855            CoinBigIndex j;
     856            for (j = columnStart[iColumn];
     857                 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     858              int iRow = row[j];
     859              value += pi[iRow] * elementByColumn[j];
     860            }
     861            if (fabs(value) > zeroTolerance) {
     862              index[numberNonZero++] = iColumn;
     863              array[iColumn] = -value;
     864            }
     865          }
     866        } else {
     867          for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
     868            double value = 0.0;
     869            CoinBigIndex j;
     870            for (j = columnStart[iColumn];
     871                 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     872              int iRow = row[j];
     873              value += pi[iRow] * elementByColumn[j];
     874            }
     875            value *= scalar;
     876            if (fabs(value) > zeroTolerance) {
     877              index[numberNonZero++] = iColumn;
     878              array[iColumn] = value;
     879            }
     880          }
     881        }
     882      } else {
    886883#ifdef CLP_INVESTIGATE
    887                     if (model->clpScaledMatrix())
    888                          printf("scaledMatrix_ at %d of ClpPackedMatrix\n", __LINE__);
    889 #endif
    890                     // scaled
    891                     if (scalar == -1.0) {
    892                          for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
    893                               double value = 0.0;
    894                               CoinBigIndex j;
    895                               const double * columnScale = model->columnScale();
    896                               for (j = columnStart[iColumn];
    897                                         j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    898                                    int iRow = row[j];
    899                                    value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
    900                               }
    901                               value *= columnScale[iColumn];
    902                               if (fabs(value) > zeroTolerance) {
    903                                    index[numberNonZero++] = iColumn;
    904                                    array[iColumn] = -value;
    905                               }
    906                          }
    907                     } else {
    908                          for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
    909                               double value = 0.0;
    910                               CoinBigIndex j;
    911                               const double * columnScale = model->columnScale();
    912                               for (j = columnStart[iColumn];
    913                                         j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    914                                    int iRow = row[j];
    915                                    value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
    916                               }
    917                               value *= scalar * columnScale[iColumn];
    918                               if (fabs(value) > zeroTolerance) {
    919                                    index[numberNonZero++] = iColumn;
    920                                    array[iColumn] = value;
    921                               }
    922                          }
    923                     }
    924                }
    925           }
    926           columnArray->setNumElements(numberNonZero);
    927           y->setNumElements(0);
    928      } else {
    929           // do by row
    930           rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray);
    931      }
    932      if (packed)
    933           columnArray->setPackedMode(true);
    934      if (0) {
    935           columnArray->checkClean();
    936           int numberNonZero = columnArray->getNumElements();;
    937           int * index = columnArray->getIndices();
    938           double * array = columnArray->denseVector();
    939           int i;
    940           for (i = 0; i < numberNonZero; i++) {
    941                int j = index[i];
    942                double value;
    943                if (packed)
    944                     value = array[i];
    945                else
    946                     value = array[j];
    947                printf("Ti %d %d %g\n", i, j, value);
    948           }
    949      }
     884        if (model->clpScaledMatrix())
     885          printf("scaledMatrix_ at %d of ClpPackedMatrix\n", __LINE__);
     886#endif
     887        // scaled
     888        if (scalar == -1.0) {
     889          for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
     890            double value = 0.0;
     891            CoinBigIndex j;
     892            const double *columnScale = model->columnScale();
     893            for (j = columnStart[iColumn];
     894                 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     895              int iRow = row[j];
     896              value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
     897            }
     898            value *= columnScale[iColumn];
     899            if (fabs(value) > zeroTolerance) {
     900              index[numberNonZero++] = iColumn;
     901              array[iColumn] = -value;
     902            }
     903          }
     904        } else {
     905          for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
     906            double value = 0.0;
     907            CoinBigIndex j;
     908            const double *columnScale = model->columnScale();
     909            for (j = columnStart[iColumn];
     910                 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     911              int iRow = row[j];
     912              value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
     913            }
     914            value *= scalar * columnScale[iColumn];
     915            if (fabs(value) > zeroTolerance) {
     916              index[numberNonZero++] = iColumn;
     917              array[iColumn] = value;
     918            }
     919          }
     920        }
     921      }
     922    }
     923    columnArray->setNumElements(numberNonZero);
     924    y->setNumElements(0);
     925  } else {
     926    // do by row
     927    rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray);
     928  }
     929  if (packed)
     930    columnArray->setPackedMode(true);
     931  if (0) {
     932    columnArray->checkClean();
     933    int numberNonZero = columnArray->getNumElements();
     934    ;
     935    int *index = columnArray->getIndices();
     936    double *array = columnArray->denseVector();
     937    int i;
     938    for (i = 0; i < numberNonZero; i++) {
     939      int j = index[i];
     940      double value;
     941      if (packed)
     942        value = array[i];
     943      else
     944        value = array[j];
     945      printf("Ti %d %d %g\n", i, j, value);
     946    }
     947  }
    950948}
    951949//static int xA=0;
     
    959957   This does by column and knows no gaps
    960958   Squashes small elements and knows about ClpSimplex */
    961 void
    962 ClpPackedMatrix::transposeTimesByColumn(const ClpSimplex * model, double scalar,
    963                                         const CoinIndexedVector * rowArray,
    964                                         CoinIndexedVector * y,
    965                                         CoinIndexedVector * columnArray) const
    966 {
    967      double * COIN_RESTRICT pi = rowArray->denseVector();
    968      int numberNonZero = 0;
    969      int * COIN_RESTRICT index = columnArray->getIndices();
    970      double * COIN_RESTRICT array = columnArray->denseVector();
    971      int numberInRowArray = rowArray->getNumElements();
    972      // maybe I need one in OsiSimplex
    973      double zeroTolerance = model->zeroTolerance();
    974      bool packed = rowArray->packedMode();
    975      // do by column
    976      int iColumn;
    977      // get matrix data pointers
    978      const int * COIN_RESTRICT row = matrix_->getIndices();
    979      const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
    980      const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
    981      const double * COIN_RESTRICT rowScale = model->rowScale();
    982      assert (!y->getNumElements());
    983      assert (numberActiveColumns_ > 0);
    984      const ClpPackedMatrix * thisMatrix = this;
     959void ClpPackedMatrix::transposeTimesByColumn(const ClpSimplex *model, double scalar,
     960  const CoinIndexedVector *rowArray,
     961  CoinIndexedVector *y,
     962  CoinIndexedVector *columnArray) const
     963{
     964  double *COIN_RESTRICT pi = rowArray->denseVector();
     965  int numberNonZero = 0;
     966  int *COIN_RESTRICT index = columnArray->getIndices();
     967  double *COIN_RESTRICT array = columnArray->denseVector();
     968  int numberInRowArray = rowArray->getNumElements();
     969  // maybe I need one in OsiSimplex
     970  double zeroTolerance = model->zeroTolerance();
     971  bool packed = rowArray->packedMode();
     972  // do by column
     973  int iColumn;
     974  // get matrix data pointers
     975  const int *COIN_RESTRICT row = matrix_->getIndices();
     976  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     977  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
     978  const double *COIN_RESTRICT rowScale = model->rowScale();
     979  assert(!y->getNumElements());
     980  assert(numberActiveColumns_ > 0);
     981  const ClpPackedMatrix *thisMatrix = this;
    985982#if 0
    986983     ClpPackedMatrix * scaledMatrix = model->clpScaledMatrix();
     
    999996     }
    1000997#endif
    1001      if (packed) {
    1002           // need to expand pi into y
    1003           assert(y->capacity() >= model->numberRows());
    1004           double * COIN_RESTRICT piOld = pi;
    1005           pi = y->denseVector();
    1006           const int * COIN_RESTRICT whichRow = rowArray->getIndices();
    1007           int i;
    1008           if (!rowScale) {
    1009                // modify pi so can collapse to one loop
    1010                if (scalar == -1.0) {
    1011                     //yA += numberInRowArray;
    1012                     for (i = 0; i < numberInRowArray; i++) {
    1013                          int iRow = whichRow[i];
    1014                          pi[iRow] = -piOld[i];
     998  if (packed) {
     999    // need to expand pi into y
     1000    assert(y->capacity() >= model->numberRows());
     1001    double *COIN_RESTRICT piOld = pi;
     1002    pi = y->denseVector();
     1003    const int *COIN_RESTRICT whichRow = rowArray->getIndices();
     1004    int i;
     1005    if (!rowScale) {
     1006      // modify pi so can collapse to one loop
     1007      if (scalar == -1.0) {
     1008        //yA += numberInRowArray;
     1009        for (i = 0; i < numberInRowArray; i++) {
     1010          int iRow = whichRow[i];
     1011          pi[iRow] = -piOld[i];
     1012        }
     1013      } else {
     1014        for (i = 0; i < numberInRowArray; i++) {
     1015          int iRow = whichRow[i];
     1016          pi[iRow] = scalar * piOld[i];
     1017        }
     1018      }
     1019      if (!columnCopy_) {
     1020        if ((model->specialOptions(), 131072) != 0) {
     1021          if (model->spareIntArray_[0] > 0) {
     1022            CoinIndexedVector *spareArray = model->rowArray(3);
     1023            // also do dualColumn stuff
     1024            double *COIN_RESTRICT spare = spareArray->denseVector();
     1025            int *COIN_RESTRICT spareIndex = spareArray->getIndices();
     1026            const double *COIN_RESTRICT reducedCost = model->djRegion(0);
     1027            double multiplier[] = { -1.0, 1.0 };
     1028            double dualT = -model->currentDualTolerance();
     1029            double acceptablePivot = model->spareDoubleArray_[0];
     1030            // We can also see if infeasible or pivoting on free
     1031            double tentativeTheta = 1.0e15;
     1032            double upperTheta = 1.0e31;
     1033            int addSequence = model->numberColumns();
     1034            const unsigned char *COIN_RESTRICT statusArray = model->statusArray() + addSequence;
     1035            int numberRemaining = 0;
     1036            assert(scalar == -1.0);
     1037            for (i = 0; i < numberInRowArray; i++) {
     1038              int iSequence = whichRow[i];
     1039              int iStatus = (statusArray[iSequence] & 3) - 1;
     1040              if (iStatus) {
     1041                double mult = multiplier[iStatus - 1];
     1042                double alpha = piOld[i] * mult;
     1043                double oldValue;
     1044                double value;
     1045                if (alpha > 0.0) {
     1046                  oldValue = reducedCost[iSequence] * mult;
     1047                  value = oldValue - tentativeTheta * alpha;
     1048                  if (value < dualT) {
     1049                    value = oldValue - upperTheta * alpha;
     1050                    if (value < dualT && alpha >= acceptablePivot) {
     1051                      upperTheta = (oldValue - dualT) / alpha;
     1052                      //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
    10151053                    }
    1016                } else {
    1017                     for (i = 0; i < numberInRowArray; i++) {
    1018                          int iRow = whichRow[i];
    1019                          pi[iRow] = scalar * piOld[i];
    1020                     }
    1021                }
    1022                if (!columnCopy_) {
    1023                     if ((model->specialOptions(), 131072) != 0) {
    1024                          if(model->spareIntArray_[0] > 0) {
    1025                               CoinIndexedVector * spareArray = model->rowArray(3);
    1026                               // also do dualColumn stuff
    1027                               double * COIN_RESTRICT spare = spareArray->denseVector();
    1028                               int * COIN_RESTRICT spareIndex = spareArray->getIndices();
    1029                               const double * COIN_RESTRICT reducedCost = model->djRegion(0);
    1030                               double multiplier[] = { -1.0, 1.0};
    1031                               double dualT = - model->currentDualTolerance();
    1032                               double acceptablePivot = model->spareDoubleArray_[0];
    1033                               // We can also see if infeasible or pivoting on free
    1034                               double tentativeTheta = 1.0e15;
    1035                               double upperTheta = 1.0e31;
    1036                               int addSequence = model->numberColumns();
    1037                               const unsigned char * COIN_RESTRICT statusArray = model->statusArray() + addSequence;
    1038                               int numberRemaining = 0;
    1039                               assert (scalar == -1.0);
    1040                               for (i = 0; i < numberInRowArray; i++) {
    1041                                    int iSequence = whichRow[i];
    1042                                    int iStatus = (statusArray[iSequence] & 3) - 1;
    1043                                    if (iStatus) {
    1044                                         double mult = multiplier[iStatus-1];
    1045                                         double alpha = piOld[i] * mult;
    1046                                         double oldValue;
    1047                                         double value;
    1048                                         if (alpha > 0.0) {
    1049                                              oldValue = reducedCost[iSequence] * mult;
    1050                                              value = oldValue - tentativeTheta * alpha;
    1051                                              if (value < dualT) {
    1052                                                   value = oldValue - upperTheta * alpha;
    1053                                                   if (value < dualT && alpha >= acceptablePivot) {
    1054                                                        upperTheta = (oldValue - dualT) / alpha;
    1055                                                        //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
    1056                                                   }
    1057                                                   // add to list
    1058                                                   spare[numberRemaining] = alpha * mult;
    1059                                                   spareIndex[numberRemaining++] = iSequence + addSequence;
    1060                                              }
    1061                                         }
    1062                                    }
    1063                               }
    1064                               numberNonZero =
    1065                                    thisMatrix->gutsOfTransposeTimesUnscaled(pi,
    1066                                              columnArray->getIndices(),
    1067                                              columnArray->denseVector(),
    1068                                              model->statusArray(),
    1069                                              spareIndex,
    1070                                              spare,
    1071                                              model->djRegion(1),
    1072                                              upperTheta,
    1073                                              acceptablePivot,
    1074                                              model->currentDualTolerance(),
    1075                                              numberRemaining,
    1076                                              zeroTolerance);
    1077                               model->spareDoubleArray_[0] = upperTheta;
    1078                               spareArray->setNumElements(numberRemaining);
     1054                    // add to list
     1055                    spare[numberRemaining] = alpha * mult;
     1056                    spareIndex[numberRemaining++] = iSequence + addSequence;
     1057                  }
     1058                }
     1059              }
     1060            }
     1061            numberNonZero = thisMatrix->gutsOfTransposeTimesUnscaled(pi,
     1062              columnArray->getIndices(),
     1063              columnArray->denseVector(),
     1064              model->statusArray(),
     1065              spareIndex,
     1066              spare,
     1067              model->djRegion(1),
     1068              upperTheta,
     1069              acceptablePivot,
     1070              model->currentDualTolerance(),
     1071              numberRemaining,
     1072              zeroTolerance);
     1073            model->spareDoubleArray_[0] = upperTheta;
     1074            spareArray->setNumElements(numberRemaining);
    10791075#if 0
    10801076                              columnArray->setNumElements(numberNonZero);
     
    10871083                                     spareArray->getNumElements(),upperTheta);
    10881084#endif
    1089                               // signal partially done
    1090                               model->spareIntArray_[0] = -2;
    1091                          } else {
    1092                               numberNonZero =
    1093                                    thisMatrix->gutsOfTransposeTimesUnscaled(pi,
    1094                                              columnArray->getIndices(),
    1095                                              columnArray->denseVector(),
    1096                                              model->statusArray(),
    1097                                              zeroTolerance);
    1098                          }
    1099                     } else {
    1100                          numberNonZero =
    1101                               thisMatrix->gutsOfTransposeTimesUnscaled(pi,
    1102                                         columnArray->getIndices(),
    1103                                         columnArray->denseVector(),
    1104                                         zeroTolerance);
    1105                     }
    1106                     columnArray->setNumElements(numberNonZero);
    1107                     //xA++;
    1108                } else {
    1109                  if ((model->moreSpecialOptions()&8)!=0 &&
    1110                      model->algorithm()<0) {
    1111                    columnCopy_->transposeTimes(model, pi, columnArray,
    1112                      model->rowArray(3),rowArray);
    1113                    model->spareIntArray_[0]=-2;
    1114                  } else {
    1115                     columnCopy_->transposeTimes(model, pi, columnArray);
    1116                  }
    1117                     numberNonZero = columnArray->getNumElements();
    1118                     //xB++;
    1119                }
     1085            // signal partially done
     1086            model->spareIntArray_[0] = -2;
    11201087          } else {
     1088            numberNonZero = thisMatrix->gutsOfTransposeTimesUnscaled(pi,
     1089              columnArray->getIndices(),
     1090              columnArray->denseVector(),
     1091              model->statusArray(),
     1092              zeroTolerance);
     1093          }
     1094        } else {
     1095          numberNonZero = thisMatrix->gutsOfTransposeTimesUnscaled(pi,
     1096            columnArray->getIndices(),
     1097            columnArray->denseVector(),
     1098            zeroTolerance);
     1099        }
     1100        columnArray->setNumElements(numberNonZero);
     1101        //xA++;
     1102      } else {
     1103        if ((model->moreSpecialOptions() & 8) != 0 && model->algorithm() < 0) {
     1104          columnCopy_->transposeTimes(model, pi, columnArray,
     1105            model->rowArray(3), rowArray);
     1106          model->spareIntArray_[0] = -2;
     1107        } else {
     1108          columnCopy_->transposeTimes(model, pi, columnArray);
     1109        }
     1110        numberNonZero = columnArray->getNumElements();
     1111        //xB++;
     1112      }
     1113    } else {
    11211114#ifdef CLP_INVESTIGATE
    1122                if (model->clpScaledMatrix())
    1123                     printf("scaledMatrix_ at %d of ClpPackedMatrix\n", __LINE__);
    1124 #endif
    1125                // scaled
    1126                // modify pi so can collapse to one loop
    1127                if (scalar == -1.0) {
    1128                     //yC += numberInRowArray;
    1129                     for (i = 0; i < numberInRowArray; i++) {
    1130                          int iRow = whichRow[i];
    1131                          pi[iRow] = -piOld[i] * rowScale[iRow];
    1132                     }
    1133                } else {
    1134                     for (i = 0; i < numberInRowArray; i++) {
    1135                          int iRow = whichRow[i];
    1136                          pi[iRow] = scalar * piOld[i] * rowScale[iRow];
    1137                     }
    1138                }
    1139                const double * columnScale = model->columnScale();
    1140                if (!columnCopy_) {
    1141                     if ((model->specialOptions(), 131072) != 0)
    1142                          numberNonZero =
    1143                               gutsOfTransposeTimesScaled(pi, columnScale,
    1144                                                          columnArray->getIndices(),
    1145                                                          columnArray->denseVector(),
    1146                                                          model->statusArray(),
    1147                                                          zeroTolerance);
    1148                     else
    1149                          numberNonZero =
    1150                               gutsOfTransposeTimesScaled(pi, columnScale,
    1151                                                          columnArray->getIndices(),
    1152                                                          columnArray->denseVector(),
    1153                                                          zeroTolerance);
    1154                     columnArray->setNumElements(numberNonZero);
    1155                     //xC++;
    1156                } else {
    1157                  if ((model->moreSpecialOptions()&8)!=0 &&
    1158                      model->algorithm()<0) {
    1159                    columnCopy_->transposeTimes(model, pi, columnArray,
    1160                      model->rowArray(3),rowArray);
    1161                    model->spareIntArray_[0]=-2;
    1162                  } else {
    1163                     columnCopy_->transposeTimes(model, pi, columnArray);
    1164                  }
    1165                     numberNonZero = columnArray->getNumElements();
    1166                     //xD++;
    1167                }
    1168           }
    1169           // zero out
    1170           int numberRows = model->numberRows();
    1171           if (numberInRowArray * 4 < numberRows) {
    1172                for (i = 0; i < numberInRowArray; i++) {
    1173                     int iRow = whichRow[i];
    1174                     pi[iRow] = 0.0;
    1175                }
    1176           } else {
    1177                CoinZeroN(pi, numberRows);
    1178           }
    1179           //int kA=xA+xB;
    1180           //int kC=xC+xD;
    1181           //if ((kA+kC)%10000==0)
    1182           //printf("AA %d %d %g, CC %d %d %g\n",
    1183           //     xA,xB,kA ? yA/(double)(kA): 0.0,xC,xD,kC ? yC/(double) (kC) :0.0);
    1184      } else {
    1185           if (!rowScale) {
    1186                if (scalar == -1.0) {
    1187                     double value = 0.0;
    1188                     CoinBigIndex j;
    1189                     CoinBigIndex end = columnStart[1];
    1190                     for (j = columnStart[0]; j < end; j++) {
    1191                          int iRow = row[j];
    1192                          value += pi[iRow] * elementByColumn[j];
    1193                     }
    1194                     for (iColumn = 0; iColumn < numberActiveColumns_ - 1; iColumn++) {
    1195                          CoinBigIndex start = end;
    1196                          end = columnStart[iColumn+2];
    1197                          if (fabs(value) > zeroTolerance) {
    1198                               array[iColumn] = -value;
    1199                               index[numberNonZero++] = iColumn;
    1200                          }
    1201                          value = 0.0;
    1202                          for (j = start; j < end; j++) {
    1203                               int iRow = row[j];
    1204                               value += pi[iRow] * elementByColumn[j];
    1205                          }
    1206                     }
    1207                     if (fabs(value) > zeroTolerance) {
    1208                          array[iColumn] = -value;
    1209                          index[numberNonZero++] = iColumn;
    1210                     }
    1211                } else {
    1212                     double value = 0.0;
    1213                     CoinBigIndex j;
    1214                     CoinBigIndex end = columnStart[1];
    1215                     for (j = columnStart[0]; j < end; j++) {
    1216                          int iRow = row[j];
    1217                          value += pi[iRow] * elementByColumn[j];
    1218                     }
    1219                     for (iColumn = 0; iColumn < numberActiveColumns_ - 1; iColumn++) {
    1220                          value *= scalar;
    1221                          CoinBigIndex start = end;
    1222                          end = columnStart[iColumn+2];
    1223                          if (fabs(value) > zeroTolerance) {
    1224                               array[iColumn] = value;
    1225                               index[numberNonZero++] = iColumn;
    1226                          }
    1227                          value = 0.0;
    1228                          for (j = start; j < end; j++) {
    1229                               int iRow = row[j];
    1230                               value += pi[iRow] * elementByColumn[j];
    1231                          }
    1232                     }
    1233                     value *= scalar;
    1234                     if (fabs(value) > zeroTolerance) {
    1235                          array[iColumn] = value;
    1236                          index[numberNonZero++] = iColumn;
    1237                     }
    1238                }
    1239           } else {
     1115      if (model->clpScaledMatrix())
     1116        printf("scaledMatrix_ at %d of ClpPackedMatrix\n", __LINE__);
     1117#endif
     1118      // scaled
     1119      // modify pi so can collapse to one loop
     1120      if (scalar == -1.0) {
     1121        //yC += numberInRowArray;
     1122        for (i = 0; i < numberInRowArray; i++) {
     1123          int iRow = whichRow[i];
     1124          pi[iRow] = -piOld[i] * rowScale[iRow];
     1125        }
     1126      } else {
     1127        for (i = 0; i < numberInRowArray; i++) {
     1128          int iRow = whichRow[i];
     1129          pi[iRow] = scalar * piOld[i] * rowScale[iRow];
     1130        }
     1131      }
     1132      const double *columnScale = model->columnScale();
     1133      if (!columnCopy_) {
     1134        if ((model->specialOptions(), 131072) != 0)
     1135          numberNonZero = gutsOfTransposeTimesScaled(pi, columnScale,
     1136            columnArray->getIndices(),
     1137            columnArray->denseVector(),
     1138            model->statusArray(),
     1139            zeroTolerance);
     1140        else
     1141          numberNonZero = gutsOfTransposeTimesScaled(pi, columnScale,
     1142            columnArray->getIndices(),
     1143            columnArray->denseVector(),
     1144            zeroTolerance);
     1145        columnArray->setNumElements(numberNonZero);
     1146        //xC++;
     1147      } else {
     1148        if ((model->moreSpecialOptions() & 8) != 0 && model->algorithm() < 0) {
     1149          columnCopy_->transposeTimes(model, pi, columnArray,
     1150            model->rowArray(3), rowArray);
     1151          model->spareIntArray_[0] = -2;
     1152        } else {
     1153          columnCopy_->transposeTimes(model, pi, columnArray);
     1154        }
     1155        numberNonZero = columnArray->getNumElements();
     1156        //xD++;
     1157      }
     1158    }
     1159    // zero out
     1160    int numberRows = model->numberRows();
     1161    if (numberInRowArray * 4 < numberRows) {
     1162      for (i = 0; i < numberInRowArray; i++) {
     1163        int iRow = whichRow[i];
     1164        pi[iRow] = 0.0;
     1165      }
     1166    } else {
     1167      CoinZeroN(pi, numberRows);
     1168    }
     1169    //int kA=xA+xB;
     1170    //int kC=xC+xD;
     1171    //if ((kA+kC)%10000==0)
     1172    //printf("AA %d %d %g, CC %d %d %g\n",
     1173    //     xA,xB,kA ? yA/(double)(kA): 0.0,xC,xD,kC ? yC/(double) (kC) :0.0);
     1174  } else {
     1175    if (!rowScale) {
     1176      if (scalar == -1.0) {
     1177        double value = 0.0;
     1178        CoinBigIndex j;
     1179        CoinBigIndex end = columnStart[1];
     1180        for (j = columnStart[0]; j < end; j++) {
     1181          int iRow = row[j];
     1182          value += pi[iRow] * elementByColumn[j];
     1183        }
     1184        for (iColumn = 0; iColumn < numberActiveColumns_ - 1; iColumn++) {
     1185          CoinBigIndex start = end;
     1186          end = columnStart[iColumn + 2];
     1187          if (fabs(value) > zeroTolerance) {
     1188            array[iColumn] = -value;
     1189            index[numberNonZero++] = iColumn;
     1190          }
     1191          value = 0.0;
     1192          for (j = start; j < end; j++) {
     1193            int iRow = row[j];
     1194            value += pi[iRow] * elementByColumn[j];
     1195          }
     1196        }
     1197        if (fabs(value) > zeroTolerance) {
     1198          array[iColumn] = -value;
     1199          index[numberNonZero++] = iColumn;
     1200        }
     1201      } else {
     1202        double value = 0.0;
     1203        CoinBigIndex j;
     1204        CoinBigIndex end = columnStart[1];
     1205        for (j = columnStart[0]; j < end; j++) {
     1206          int iRow = row[j];
     1207          value += pi[iRow] * elementByColumn[j];
     1208        }
     1209        for (iColumn = 0; iColumn < numberActiveColumns_ - 1; iColumn++) {
     1210          value *= scalar;
     1211          CoinBigIndex start = end;
     1212          end = columnStart[iColumn + 2];
     1213          if (fabs(value) > zeroTolerance) {
     1214            array[iColumn] = value;
     1215            index[numberNonZero++] = iColumn;
     1216          }
     1217          value = 0.0;
     1218          for (j = start; j < end; j++) {
     1219            int iRow = row[j];
     1220            value += pi[iRow] * elementByColumn[j];
     1221          }
     1222        }
     1223        value *= scalar;
     1224        if (fabs(value) > zeroTolerance) {
     1225          array[iColumn] = value;
     1226          index[numberNonZero++] = iColumn;
     1227        }
     1228      }
     1229    } else {
    12401230#ifdef CLP_INVESTIGATE
    1241                if (model->clpScaledMatrix())
    1242                     printf("scaledMatrix_ at %d of ClpPackedMatrix\n", __LINE__);
    1243 #endif
    1244                // scaled
    1245                if (scalar == -1.0) {
    1246                     const double * COIN_RESTRICT columnScale = model->columnScale();
    1247                     double value = 0.0;
    1248                     double scale = columnScale[0];
    1249                     CoinBigIndex j;
    1250                     CoinBigIndex end = columnStart[1];
    1251                     for (j = columnStart[0]; j < end; j++) {
    1252                          int iRow = row[j];
    1253                          value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
    1254                     }
    1255                     for (iColumn = 0; iColumn < numberActiveColumns_ - 1; iColumn++) {
    1256                          value *= scale;
    1257                          CoinBigIndex start = end;
    1258                          end = columnStart[iColumn+2];
    1259                          scale = columnScale[iColumn+1];
    1260                          if (fabs(value) > zeroTolerance) {
    1261                               array[iColumn] = -value;
    1262                               index[numberNonZero++] = iColumn;
    1263                          }
    1264                          value = 0.0;
    1265                          for (j = start; j < end; j++) {
    1266                               int iRow = row[j];
    1267                               value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
    1268                          }
    1269                     }
    1270                     value *= scale;
    1271                     if (fabs(value) > zeroTolerance) {
    1272                          array[iColumn] = -value;
    1273                          index[numberNonZero++] = iColumn;
    1274                     }
    1275                } else {
    1276                     const double * COIN_RESTRICT columnScale = model->columnScale();
    1277                     double value = 0.0;
    1278                     double scale = columnScale[0] * scalar;
    1279                     CoinBigIndex j;
    1280                     CoinBigIndex end = columnStart[1];
    1281                     for (j = columnStart[0]; j < end; j++) {
    1282                          int iRow = row[j];
    1283                          value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
    1284                     }
    1285                     for (iColumn = 0; iColumn < numberActiveColumns_ - 1; iColumn++) {
    1286                          value *= scale;
    1287                          CoinBigIndex start = end;
    1288                          end = columnStart[iColumn+2];
    1289                          scale = columnScale[iColumn+1] * scalar;
    1290                          if (fabs(value) > zeroTolerance) {
    1291                               array[iColumn] = value;
    1292                               index[numberNonZero++] = iColumn;
    1293                          }
    1294                          value = 0.0;
    1295                          for (j = start; j < end; j++) {
    1296                               int iRow = row[j];
    1297                               value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
    1298                          }
    1299                     }
    1300                     value *= scale;
    1301                     if (fabs(value) > zeroTolerance) {
    1302                          array[iColumn] = value;
    1303                          index[numberNonZero++] = iColumn;
    1304                     }
    1305                }
    1306           }
    1307      }
    1308      columnArray->setNumElements(numberNonZero);
    1309      y->setNumElements(0);
    1310      if (packed)
    1311           columnArray->setPackedMode(true);
     1231      if (model->clpScaledMatrix())
     1232        printf("scaledMatrix_ at %d of ClpPackedMatrix\n", __LINE__);
     1233#endif
     1234      // scaled
     1235      if (scalar == -1.0) {
     1236        const double *COIN_RESTRICT columnScale = model->columnScale();
     1237        double value = 0.0;
     1238        double scale = columnScale[0];
     1239        CoinBigIndex j;
     1240        CoinBigIndex end = columnStart[1];
     1241        for (j = columnStart[0]; j < end; j++) {
     1242          int iRow = row[j];
     1243          value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
     1244        }
     1245        for (iColumn = 0; iColumn < numberActiveColumns_ - 1; iColumn++) {
     1246          value *= scale;
     1247          CoinBigIndex start = end;
     1248          end = columnStart[iColumn + 2];
     1249          scale = columnScale[iColumn + 1];
     1250          if (fabs(value) > zeroTolerance) {
     1251            array[iColumn] = -value;
     1252            index[numberNonZero++] = iColumn;
     1253          }
     1254          value = 0.0;
     1255          for (j = start; j < end; j++) {
     1256            int iRow = row[j];
     1257            value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
     1258          }
     1259        }
     1260        value *= scale;
     1261        if (fabs(value) > zeroTolerance) {
     1262          array[iColumn] = -value;
     1263          index[numberNonZero++] = iColumn;
     1264        }
     1265      } else {
     1266        const double *COIN_RESTRICT columnScale = model->columnScale();
     1267        double value = 0.0;
     1268        double scale = columnScale[0] * scalar;
     1269        CoinBigIndex j;
     1270        CoinBigIndex end = columnStart[1];
     1271        for (j = columnStart[0]; j < end; j++) {
     1272          int iRow = row[j];
     1273          value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
     1274        }
     1275        for (iColumn = 0; iColumn < numberActiveColumns_ - 1; iColumn++) {
     1276          value *= scale;
     1277          CoinBigIndex start = end;
     1278          end = columnStart[iColumn + 2];
     1279          scale = columnScale[iColumn + 1] * scalar;
     1280          if (fabs(value) > zeroTolerance) {
     1281            array[iColumn] = value;
     1282            index[numberNonZero++] = iColumn;
     1283          }
     1284          value = 0.0;
     1285          for (j = start; j < end; j++) {
     1286            int iRow = row[j];
     1287            value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
     1288          }
     1289        }
     1290        value *= scale;
     1291        if (fabs(value) > zeroTolerance) {
     1292          array[iColumn] = value;
     1293          index[numberNonZero++] = iColumn;
     1294        }
     1295      }
     1296    }
     1297  }
     1298  columnArray->setNumElements(numberNonZero);
     1299  y->setNumElements(0);
     1300  if (packed)
     1301    columnArray->setPackedMode(true);
    13121302}
    13131303/* Return <code>x * A in <code>z</code>.
    13141304        Squashes small elements and knows about ClpSimplex */
    1315 void
    1316 ClpPackedMatrix::transposeTimesByRow(const ClpSimplex * model, double scalar,
    1317                                      const CoinIndexedVector * rowArray,
    1318                                      CoinIndexedVector * y,
    1319                                      CoinIndexedVector * columnArray) const
    1320 {
    1321      columnArray->clear();
    1322      double * COIN_RESTRICT pi = rowArray->denseVector();
    1323      int numberNonZero = 0;
    1324      int * COIN_RESTRICT index = columnArray->getIndices();
    1325      double * COIN_RESTRICT array = columnArray->denseVector();
    1326      int numberInRowArray = rowArray->getNumElements();
    1327      // maybe I need one in OsiSimplex
    1328      double zeroTolerance = model->zeroTolerance();
    1329      const int * COIN_RESTRICT column = matrix_->getIndices();
    1330      const CoinBigIndex * COIN_RESTRICT rowStart = getVectorStarts();
    1331      const double * COIN_RESTRICT element = getElements();
    1332      const int * COIN_RESTRICT whichRow = rowArray->getIndices();
    1333      bool packed = rowArray->packedMode();
    1334      if (numberInRowArray > 2) {
    1335           // do by rows
    1336           // ** Row copy is already scaled
    1337           int iRow;
    1338           int i;
    1339           int numberOriginal = 0;
    1340           if (packed) {
    1341                int * COIN_RESTRICT index = columnArray->getIndices();
    1342                double * COIN_RESTRICT array = columnArray->denseVector();
     1305void ClpPackedMatrix::transposeTimesByRow(const ClpSimplex *model, double scalar,
     1306  const CoinIndexedVector *rowArray,
     1307  CoinIndexedVector *y,
     1308  CoinIndexedVector *columnArray) const
     1309{
     1310  columnArray->clear();
     1311  double *COIN_RESTRICT pi = rowArray->denseVector();
     1312  int numberNonZero = 0;
     1313  int *COIN_RESTRICT index = columnArray->getIndices();
     1314  double *COIN_RESTRICT array = columnArray->denseVector();
     1315  int numberInRowArray = rowArray->getNumElements();
     1316  // maybe I need one in OsiSimplex
     1317  double zeroTolerance = model->zeroTolerance();
     1318  const int *COIN_RESTRICT column = matrix_->getIndices();
     1319  const CoinBigIndex *COIN_RESTRICT rowStart = getVectorStarts();
     1320  const double *COIN_RESTRICT element = getElements();
     1321  const int *COIN_RESTRICT whichRow = rowArray->getIndices();
     1322  bool packed = rowArray->packedMode();
     1323  if (numberInRowArray > 2) {
     1324    // do by rows
     1325    // ** Row copy is already scaled
     1326    int iRow;
     1327    int i;
     1328    int numberOriginal = 0;
     1329    if (packed) {
     1330      int *COIN_RESTRICT index = columnArray->getIndices();
     1331      double *COIN_RESTRICT array = columnArray->denseVector();
    13431332#if 0
    13441333               {
     
    13521341#endif
    13531342#define COIN_SPARSE_MATRIX 2
    1354                CoinBigIndex numberCovered=0;
    1355                int numberColumns = matrix_->getNumCols();
    1356                bool sparse=true;
    1357                for (int i = 0; i < numberInRowArray; i++) {
    1358                 int iRow = whichRow[i];
    1359                  numberCovered += rowStart[iRow+1] - rowStart[iRow];
    1360                 // ? exact ratio
    1361                  if (numberCovered>numberColumns) {
    1362                    sparse=false;
    1363                    break;
    1364                 }
    1365                }
    1366                if (sparse) {
    1367                  assert (!y->getNumElements());
     1343      CoinBigIndex numberCovered = 0;
     1344      int numberColumns = matrix_->getNumCols();
     1345      bool sparse = true;
     1346      for (int i = 0; i < numberInRowArray; i++) {
     1347        int iRow = whichRow[i];
     1348        numberCovered += rowStart[iRow + 1] - rowStart[iRow];
     1349        // ? exact ratio
     1350        if (numberCovered > numberColumns) {
     1351          sparse = false;
     1352          break;
     1353        }
     1354      }
     1355      if (sparse) {
     1356        assert(!y->getNumElements());
    13681357#if COIN_SPARSE_MATRIX != 2
    1369                 // and set up mark as char array
    1370                  char * COIN_RESTRICT marked = reinterpret_cast<char *> (index+columnArray->capacity());
    1371                  int * COIN_RESTRICT lookup = y->getIndices();
     1358        // and set up mark as char array
     1359        char *COIN_RESTRICT marked = reinterpret_cast< char * >(index + columnArray->capacity());
     1360        int *COIN_RESTRICT lookup = y->getIndices();
    13721361#ifndef NDEBUG
    1373                 //int numberColumns = matrix_->getNumCols();
    1374                 //for (int i=0;i<numberColumns;i++)
    1375                 //assert(!marked[i]);
    1376 #endif
    1377                  numberNonZero=gutsOfTransposeTimesByRowGE3a(rowArray,index,array,
    1378                                                              lookup,marked,zeroTolerance,scalar);
     1362        //int numberColumns = matrix_->getNumCols();
     1363        //for (int i=0;i<numberColumns;i++)
     1364        //assert(!marked[i]);
     1365#endif
     1366        numberNonZero = gutsOfTransposeTimesByRowGE3a(rowArray, index, array,
     1367          lookup, marked, zeroTolerance, scalar);
    13791368#else
    1380                  double  * COIN_RESTRICT array2 = y->denseVector();
    1381                  numberNonZero=gutsOfTransposeTimesByRowGE3(rowArray,index,array,
    1382                                                             array2,zeroTolerance,scalar);
    1383 #endif
    1384                } else {
    1385                 numberNonZero = gutsOfTransposeTimesByRowGEK(rowArray, index, array,
    1386                                                               numberColumns, zeroTolerance, scalar);
    1387                }
    1388                columnArray->setNumElements(numberNonZero);
    1389           } else {
    1390                double * COIN_RESTRICT markVector = y->denseVector();
    1391                numberNonZero = 0;
    1392                // and set up mark as char array
    1393                char * COIN_RESTRICT marked = reinterpret_cast<char *> (markVector);
    1394                for (i = 0; i < numberOriginal; i++) {
    1395                     int iColumn = index[i];
    1396                     marked[iColumn] = 0;
    1397                }
     1369        double *COIN_RESTRICT array2 = y->denseVector();
     1370        numberNonZero = gutsOfTransposeTimesByRowGE3(rowArray, index, array,
     1371          array2, zeroTolerance, scalar);
     1372#endif
     1373      } else {
     1374        numberNonZero = gutsOfTransposeTimesByRowGEK(rowArray, index, array,
     1375          numberColumns, zeroTolerance, scalar);
     1376      }
     1377      columnArray->setNumElements(numberNonZero);
     1378    } else {
     1379      double *COIN_RESTRICT markVector = y->denseVector();
     1380      numberNonZero = 0;
     1381      // and set up mark as char array
     1382      char *COIN_RESTRICT marked = reinterpret_cast< char * >(markVector);
     1383      for (i = 0; i < numberOriginal; i++) {
     1384        int iColumn = index[i];
     1385        marked[iColumn] = 0;
     1386      }
    13981387
    1399                for (i = 0; i < numberInRowArray; i++) {
    1400                     iRow = whichRow[i];
    1401                     double value = pi[iRow] * scalar;
    1402                     CoinBigIndex j;
    1403                     for (j = rowStart[iRow]; j < rowStart[iRow+1]; j++) {
    1404                          int iColumn = column[j];
    1405                          if (!marked[iColumn]) {
    1406                               marked[iColumn] = 1;
    1407                               index[numberNonZero++] = iColumn;
    1408                          }
    1409                          array[iColumn] += value * element[j];
    1410                     }
    1411                }
    1412                // get rid of tiny values and zero out marked
    1413                numberOriginal = numberNonZero;
    1414                numberNonZero = 0;
    1415                for (i = 0; i < numberOriginal; i++) {
    1416                     int iColumn = index[i];
    1417                     marked[iColumn] = 0;
    1418                     if (fabs(array[iColumn]) > zeroTolerance) {
    1419                          index[numberNonZero++] = iColumn;
    1420                     } else {
    1421                          array[iColumn] = 0.0;
    1422                     }
    1423                }
    1424           }
    1425      } else if (numberInRowArray == 2) {
    1426           // do by rows when two rows
    1427           int numberOriginal;
    1428           int i;
    1429           CoinBigIndex j;
    1430           numberNonZero = 0;
     1388      for (i = 0; i < numberInRowArray; i++) {
     1389        iRow = whichRow[i];
     1390        double value = pi[iRow] * scalar;
     1391        CoinBigIndex j;
     1392        for (j = rowStart[iRow]; j < rowStart[iRow + 1]; j++) {
     1393          int iColumn = column[j];
     1394          if (!marked[iColumn]) {
     1395            marked[iColumn] = 1;
     1396            index[numberNonZero++] = iColumn;
     1397          }
     1398          array[iColumn] += value * element[j];
     1399        }
     1400      }
     1401      // get rid of tiny values and zero out marked
     1402      numberOriginal = numberNonZero;
     1403      numberNonZero = 0;
     1404      for (i = 0; i < numberOriginal; i++) {
     1405        int iColumn = index[i];
     1406        marked[iColumn] = 0;
     1407        if (fabs(array[iColumn]) > zeroTolerance) {
     1408          index[numberNonZero++] = iColumn;
     1409        } else {
     1410          array[iColumn] = 0.0;
     1411        }
     1412      }
     1413    }
     1414  } else if (numberInRowArray == 2) {
     1415    // do by rows when two rows
     1416    int numberOriginal;
     1417    int i;
     1418    CoinBigIndex j;
     1419    numberNonZero = 0;
    14311420
    1432           double value;
    1433           if (packed) {
    1434                gutsOfTransposeTimesByRowEQ2(rowArray, columnArray, y, zeroTolerance, scalar);
    1435                numberNonZero = columnArray->getNumElements();
    1436           } else {
    1437                int iRow = whichRow[0];
    1438                value = pi[iRow] * scalar;
    1439                for (j = rowStart[iRow]; j < rowStart[iRow+1]; j++) {
    1440                     int iColumn = column[j];
    1441                     double value2 = value * element[j];
    1442                     index[numberNonZero++] = iColumn;
    1443                     array[iColumn] = value2;
    1444                }
    1445                iRow = whichRow[1];
    1446                value = pi[iRow] * scalar;
    1447                for (j = rowStart[iRow]; j < rowStart[iRow+1]; j++) {
    1448                     int iColumn = column[j];
    1449                     double value2 = value * element[j];
    1450                     // I am assuming no zeros in matrix
    1451                     if (array[iColumn])
    1452                          value2 += array[iColumn];
    1453                     else
    1454                          index[numberNonZero++] = iColumn;
    1455                     array[iColumn] = value2;
    1456                }
    1457                // get rid of tiny values and zero out marked
    1458                numberOriginal = numberNonZero;
    1459                numberNonZero = 0;
    1460                for (i = 0; i < numberOriginal; i++) {
    1461                     int iColumn = index[i];
    1462                     if (fabs(array[iColumn]) > zeroTolerance) {
    1463                          index[numberNonZero++] = iColumn;
    1464                     } else {
    1465                          array[iColumn] = 0.0;
    1466                     }
    1467                }
    1468           }
    1469      } else if (numberInRowArray == 1) {
    1470           // Just one row
    1471           int iRow = rowArray->getIndices()[0];
    1472           numberNonZero = 0;
    1473           CoinBigIndex j;
    1474           if (packed) {
    1475                gutsOfTransposeTimesByRowEQ1(rowArray, columnArray, zeroTolerance, scalar);
    1476                numberNonZero = columnArray->getNumElements();
    1477           } else {
    1478                double value = pi[iRow] * scalar;
    1479                for (j = rowStart[iRow]; j < rowStart[iRow+1]; j++) {
    1480                     int iColumn = column[j];
    1481                     double value2 = value * element[j];
    1482                     if (fabs(value2) > zeroTolerance) {
    1483                          index[numberNonZero++] = iColumn;
    1484                          array[iColumn] = value2;
    1485                     }
    1486                }
    1487           }
    1488      }
    1489      columnArray->setNumElements(numberNonZero);
    1490      y->setNumElements(0);
     1421    double value;
     1422    if (packed) {
     1423      gutsOfTransposeTimesByRowEQ2(rowArray, columnArray, y, zeroTolerance, scalar);
     1424      numberNonZero = columnArray->getNumElements();
     1425    } else {
     1426      int iRow = whichRow[0];
     1427      value = pi[iRow] * scalar;
     1428      for (j = rowStart[iRow]; j < rowStart[iRow + 1]; j++) {
     1429        int iColumn = column[j];
     1430        double value2 = value * element[j];
     1431        index[numberNonZero++] = iColumn;
     1432        array[iColumn] = value2;
     1433      }
     1434      iRow = whichRow[1];
     1435      value = pi[iRow] * scalar;
     1436      for (j = rowStart[iRow]; j < rowStart[iRow + 1]; j++) {
     1437        int iColumn = column[j];
     1438        double value2 = value * element[j];
     1439        // I am assuming no zeros in matrix
     1440        if (array[iColumn])
     1441          value2 += array[iColumn];
     1442        else
     1443          index[numberNonZero++] = iColumn;
     1444        array[iColumn] = value2;
     1445      }
     1446      // get rid of tiny values and zero out marked
     1447      numberOriginal = numberNonZero;
     1448      numberNonZero = 0;
     1449      for (i = 0; i < numberOriginal; i++) {
     1450        int iColumn = index[i];
     1451        if (fabs(array[iColumn]) > zeroTolerance) {
     1452          index[numberNonZero++] = iColumn;
     1453        } else {
     1454          array[iColumn] = 0.0;
     1455        }
     1456      }
     1457    }
     1458  } else if (numberInRowArray == 1) {
     1459    // Just one row
     1460    int iRow = rowArray->getIndices()[0];
     1461    numberNonZero = 0;
     1462    CoinBigIndex j;
     1463    if (packed) {
     1464      gutsOfTransposeTimesByRowEQ1(rowArray, columnArray, zeroTolerance, scalar);
     1465      numberNonZero = columnArray->getNumElements();
     1466    } else {
     1467      double value = pi[iRow] * scalar;
     1468      for (j = rowStart[iRow]; j < rowStart[iRow + 1]; j++) {
     1469        int iColumn = column[j];
     1470        double value2 = value * element[j];
     1471        if (fabs(value2) > zeroTolerance) {
     1472          index[numberNonZero++] = iColumn;
     1473          array[iColumn] = value2;
     1474        }
     1475      }
     1476    }
     1477  }
     1478  columnArray->setNumElements(numberNonZero);
     1479  y->setNumElements(0);
    14911480}
    14921481// Meat of transposeTimes by column when not scaled
    1493 int
    1494 ClpPackedMatrix::gutsOfTransposeTimesUnscaled(const double * COIN_RESTRICT pi,
    1495           int * COIN_RESTRICT index,
    1496           double * COIN_RESTRICT array,
    1497           const double zeroTolerance) const
    1498 {
    1499      int numberNonZero = 0;
    1500      // get matrix data pointers
    1501      const int * COIN_RESTRICT row = matrix_->getIndices();
    1502      const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
    1503      const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
     1482int ClpPackedMatrix::gutsOfTransposeTimesUnscaled(const double *COIN_RESTRICT pi,
     1483  int *COIN_RESTRICT index,
     1484  double *COIN_RESTRICT array,
     1485  const double zeroTolerance) const
     1486{
     1487  int numberNonZero = 0;
     1488  // get matrix data pointers
     1489  const int *COIN_RESTRICT row = matrix_->getIndices();
     1490  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     1491  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
    15041492#if 1 //ndef INTEL_MKL
    1505      double value = 0.0;
    1506      CoinBigIndex j;
    1507      CoinBigIndex end = columnStart[1];
    1508      for (j = columnStart[0]; j < end; j++) {
    1509           int iRow = row[j];
    1510           value += pi[iRow] * elementByColumn[j];
    1511      }
    1512      int iColumn;
    1513      for (iColumn = 0; iColumn < numberActiveColumns_ - 1; iColumn++) {
    1514           CoinBigIndex start = end;
    1515           end = columnStart[iColumn+2];
    1516           if (fabs(value) > zeroTolerance) {
    1517                array[numberNonZero] = value;
    1518                index[numberNonZero++] = iColumn;
    1519           }
    1520           value = 0.0;
    1521           for (j = start; j < end; j++) {
    1522                int iRow = row[j];
    1523                value += pi[iRow] * elementByColumn[j];
    1524           }
    1525      }
    1526      if (fabs(value) > zeroTolerance) {
    1527           array[numberNonZero] = value;
    1528           index[numberNonZero++] = iColumn;
    1529      }
     1493  double value = 0.0;
     1494  CoinBigIndex j;
     1495  CoinBigIndex end = columnStart[1];
     1496  for (j = columnStart[0]; j < end; j++) {
     1497    int iRow = row[j];
     1498    value += pi[iRow] * elementByColumn[j];
     1499  }
     1500  int iColumn;
     1501  for (iColumn = 0; iColumn < numberActiveColumns_ - 1; iColumn++) {
     1502    CoinBigIndex start = end;
     1503    end = columnStart[iColumn + 2];
     1504    if (fabs(value) > zeroTolerance) {
     1505      array[numberNonZero] = value;
     1506      index[numberNonZero++] = iColumn;
     1507    }
     1508    value = 0.0;
     1509    for (j = start; j < end; j++) {
     1510      int iRow = row[j];
     1511      value += pi[iRow] * elementByColumn[j];
     1512    }
     1513  }
     1514  if (fabs(value) > zeroTolerance) {
     1515    array[numberNonZero] = value;
     1516    index[numberNonZero++] = iColumn;
     1517  }
    15301518#else
    1531      char transA = 'N';
    1532      //int numberRows = matrix_->getNumRows();
    1533      mkl_cspblas_dcsrgemv(&transA, const_cast<int *>(&numberActiveColumns_),
    1534                           const_cast<double *>(elementByColumn),
    1535                           const_cast<int *>(columnStart),
    1536                           const_cast<int *>(row),
    1537                           const_cast<double *>(pi), array);
    1538      int iColumn;
    1539      for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
    1540           double value = array[iColumn];
    1541           if (value) {
    1542                array[iColumn] = 0.0;
    1543                if (fabs(value) > zeroTolerance) {
    1544                     array[numberNonZero] = value;
    1545                     index[numberNonZero++] = iColumn;
    1546                }
    1547           }
    1548      }
    1549 #endif
    1550      return numberNonZero;
     1519  char transA = 'N';
     1520  //int numberRows = matrix_->getNumRows();
     1521  mkl_cspblas_dcsrgemv(&transA, const_cast< int * >(&numberActiveColumns_),
     1522    const_cast< double * >(elementByColumn),
     1523    const_cast< int * >(columnStart),
     1524    const_cast< int * >(row),
     1525    const_cast< double * >(pi), array);
     1526  int iColumn;
     1527  for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
     1528    double value = array[iColumn];
     1529    if (value) {
     1530      array[iColumn] = 0.0;
     1531      if (fabs(value) > zeroTolerance) {
     1532        array[numberNonZero] = value;
     1533        index[numberNonZero++] = iColumn;
     1534      }
     1535    }
     1536  }
     1537#endif
     1538  return numberNonZero;
    15511539}
    15521540// Meat of transposeTimes by column when scaled
    1553 int
    1554 ClpPackedMatrix::gutsOfTransposeTimesScaled(const double * COIN_RESTRICT pi,
    1555           const double * COIN_RESTRICT columnScale,
    1556           int * COIN_RESTRICT index,
    1557           double * COIN_RESTRICT array,
    1558           const double zeroTolerance) const
    1559 {
    1560      int numberNonZero = 0;
    1561      // get matrix data pointers
    1562      const int * COIN_RESTRICT row = matrix_->getIndices();
    1563      const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
    1564      const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
    1565      double value = 0.0;
    1566      double scale = columnScale[0];
    1567      CoinBigIndex j;
    1568      CoinBigIndex end = columnStart[1];
    1569      for (j = columnStart[0]; j < end; j++) {
    1570           int iRow = row[j];
    1571           value += pi[iRow] * elementByColumn[j];
    1572      }
    1573      int iColumn;
    1574      for (iColumn = 0; iColumn < numberActiveColumns_ - 1; iColumn++) {
    1575           value *= scale;
    1576           CoinBigIndex start = end;
    1577           scale = columnScale[iColumn+1];
    1578           end = columnStart[iColumn+2];
    1579           if (fabs(value) > zeroTolerance) {
    1580                array[numberNonZero] = value;
    1581                index[numberNonZero++] = iColumn;
    1582           }
    1583           value = 0.0;
    1584           for (j = start; j < end; j++) {
    1585                int iRow = row[j];
    1586                value += pi[iRow] * elementByColumn[j];
    1587           }
    1588      }
    1589      value *= scale;
    1590      if (fabs(value) > zeroTolerance) {
    1591           array[numberNonZero] = value;
    1592           index[numberNonZero++] = iColumn;
    1593      }
    1594      return numberNonZero;
     1541int ClpPackedMatrix::gutsOfTransposeTimesScaled(const double *COIN_RESTRICT pi,
     1542  const double *COIN_RESTRICT columnScale,
     1543  int *COIN_RESTRICT index,
     1544  double *COIN_RESTRICT array,
     1545  const double zeroTolerance) const
     1546{
     1547  int numberNonZero = 0;
     1548  // get matrix data pointers
     1549  const int *COIN_RESTRICT row = matrix_->getIndices();
     1550  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     1551  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
     1552  double value = 0.0;
     1553  double scale = columnScale[0];
     1554  CoinBigIndex j;
     1555  CoinBigIndex end = columnStart[1];
     1556  for (j = columnStart[0]; j < end; j++) {
     1557    int iRow = row[j];
     1558    value += pi[iRow] * elementByColumn[j];
     1559  }
     1560  int iColumn;
     1561  for (iColumn = 0; iColumn < numberActiveColumns_ - 1; iColumn++) {
     1562    value *= scale;
     1563    CoinBigIndex start = end;
     1564    scale = columnScale[iColumn + 1];
     1565    end = columnStart[iColumn + 2];
     1566    if (fabs(value) > zeroTolerance) {
     1567      array[numberNonZero] = value;
     1568      index[numberNonZero++] = iColumn;
     1569    }
     1570    value = 0.0;
     1571    for (j = start; j < end; j++) {
     1572      int iRow = row[j];
     1573      value += pi[iRow] * elementByColumn[j];
     1574    }
     1575  }
     1576  value *= scale;
     1577  if (fabs(value) > zeroTolerance) {
     1578    array[numberNonZero] = value;
     1579    index[numberNonZero++] = iColumn;
     1580  }
     1581  return numberNonZero;
    15951582}
    15961583#if ABOCA_LITE
    15971584static void
    1598 transposeTimesUnscaledBit(clpTempInfo & info)
    1599 {
    1600   const CoinBigIndex * COIN_RESTRICT columnStart=info.start;
    1601   const int * COIN_RESTRICT row = info.row;
    1602   const double * COIN_RESTRICT elementByColumn = info.element;
    1603   double * COIN_RESTRICT array = info.infeas;
    1604   int * COIN_RESTRICT index = info.which;
    1605   const unsigned char * COIN_RESTRICT status = info.status;
    1606   double zeroTolerance=info.tolerance;
    1607   const double * COIN_RESTRICT pi = info.work;
     1585transposeTimesUnscaledBit(clpTempInfo &info)
     1586{
     1587  const CoinBigIndex *COIN_RESTRICT columnStart = info.start;
     1588  const int *COIN_RESTRICT row = info.row;
     1589  const double *COIN_RESTRICT elementByColumn = info.element;
     1590  double *COIN_RESTRICT array = info.infeas;
     1591  int *COIN_RESTRICT index = info.which;
     1592  const unsigned char *COIN_RESTRICT status = info.status;
     1593  double zeroTolerance = info.tolerance;
     1594  const double *COIN_RESTRICT pi = info.work;
    16081595  int first = info.startColumn;
    1609   int last = first+info.numberToDo;
     1596  int last = first + info.numberToDo;
    16101597  double value = 0.0;
    16111598  int jColumn = -1;
    1612   int numberNonZero=0;
     1599  int numberNonZero = 0;
    16131600  for (int iColumn = first; iColumn < last; iColumn++) {
    16141601    bool wanted = ((status[iColumn] & 3) != 1);
     
    16201607    if (wanted) {
    16211608      CoinBigIndex start = columnStart[iColumn];
    1622       CoinBigIndex end = columnStart[iColumn+1];
     1609      CoinBigIndex end = columnStart[iColumn + 1];
    16231610      jColumn = iColumn;
    1624       int n = static_cast<int>(end - start);
     1611      int n = static_cast< int >(end - start);
    16251612      bool odd = (n & 1) != 0;
    16261613      n = n >> 1;
    1627       const int * COIN_RESTRICT rowThis = row + start;
    1628       const double * COIN_RESTRICT elementThis = elementByColumn + start;
     1614      const int *COIN_RESTRICT rowThis = row + start;
     1615      const double *COIN_RESTRICT elementThis = elementByColumn + start;
    16291616      for (; n; n--) {
    1630         int iRow0 = *rowThis;
    1631         int iRow1 = *(rowThis + 1);
    1632         rowThis += 2;
    1633         value += pi[iRow0] * (*elementThis);
    1634         value += pi[iRow1] * (*(elementThis + 1));
    1635         elementThis += 2;
     1617        int iRow0 = *rowThis;
     1618        int iRow1 = *(rowThis + 1);
     1619        rowThis += 2;
     1620        value += pi[iRow0] * (*elementThis);
     1621        value += pi[iRow1] * (*(elementThis + 1));
     1622        elementThis += 2;
    16361623      }
    16371624      if (odd) {
    1638         int iRow = *rowThis;
    1639         value += pi[iRow] * (*elementThis);
     1625        int iRow = *rowThis;
     1626        value += pi[iRow] * (*elementThis);
    16401627      }
    16411628    }
     
    16451632    index[numberNonZero++] = jColumn;
    16461633  }
    1647   info.numberAdded=numberNonZero;
     1634  info.numberAdded = numberNonZero;
    16481635}
    16491636#endif
    16501637// Meat of transposeTimes by column when not scaled
    1651 int
    1652 ClpPackedMatrix::gutsOfTransposeTimesUnscaled(const double * COIN_RESTRICT pi,
    1653           int * COIN_RESTRICT index,
    1654           double * COIN_RESTRICT array,
    1655           const unsigned char * COIN_RESTRICT status,
    1656           const double zeroTolerance) const
    1657 {
    1658      int numberNonZero = 0;
    1659      // get matrix data pointers
    1660      const int * COIN_RESTRICT row = matrix_->getIndices();
    1661      const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
    1662      const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
     1638int ClpPackedMatrix::gutsOfTransposeTimesUnscaled(const double *COIN_RESTRICT pi,
     1639  int *COIN_RESTRICT index,
     1640  double *COIN_RESTRICT array,
     1641  const unsigned char *COIN_RESTRICT status,
     1642  const double zeroTolerance) const
     1643{
     1644  int numberNonZero = 0;
     1645  // get matrix data pointers
     1646  const int *COIN_RESTRICT row = matrix_->getIndices();
     1647  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     1648  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
    16631649#if ABOCA_LITE
    1664      int numberThreads=abcState();
    1665      if (numberThreads) {
    1666      clpTempInfo info[ABOCA_LITE];
    1667      int chunk = (numberActiveColumns_+numberThreads-1)/numberThreads;
    1668      int n=0;
    1669      for (int i=0;i<numberThreads;i++) {
    1670        info[i].which = index+n;
    1671        info[i].infeas = array+n;
    1672        info[i].work=const_cast<double *>(pi);
    1673        info[i].status=const_cast<unsigned char *>(status);
    1674        info[i].startColumn=n;
    1675        info[i].element=elementByColumn;
    1676        info[i].start=columnStart;
    1677        info[i].row=row;
    1678        info[i].numberToDo=CoinMin(chunk,numberActiveColumns_-n);
    1679        info[i].tolerance=zeroTolerance;
    1680        n += chunk;
    1681      }
    1682      for (int i=0;i<numberThreads;i++) {
    1683        cilk_spawn transposeTimesUnscaledBit(info[i]);
    1684      }
    1685      cilk_sync;
    1686      for (int i=0;i<numberThreads;i++)
    1687        numberNonZero += info[i].numberAdded;
    1688      moveAndZero(info,2,NULL);
    1689        } else {
    1690 #endif
    1691      double value = 0.0;
    1692      int jColumn = -1;
    1693      for (int iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
    1694           bool wanted = ((status[iColumn] & 3) != 1);
    1695           if (fabs(value) > zeroTolerance) {
    1696                array[numberNonZero] = value;
    1697                index[numberNonZero++] = jColumn;
    1698           }
    1699           value = 0.0;
    1700           if (wanted) {
    1701                CoinBigIndex start = columnStart[iColumn];
    1702                CoinBigIndex end = columnStart[iColumn+1];
    1703                jColumn = iColumn;
    1704                int n = static_cast<int>(end - start);
    1705                bool odd = (n & 1) != 0;
    1706                n = n >> 1;
    1707                const int * COIN_RESTRICT rowThis = row + start;
    1708                const double * COIN_RESTRICT elementThis = elementByColumn + start;
    1709                for (; n; n--) {
    1710                     int iRow0 = *rowThis;
    1711                     int iRow1 = *(rowThis + 1);
    1712                     rowThis += 2;
    1713                     value += pi[iRow0] * (*elementThis);
    1714                     value += pi[iRow1] * (*(elementThis + 1));
    1715                     elementThis += 2;
    1716                }
    1717                if (odd) {
    1718                     int iRow = *rowThis;
    1719                     value += pi[iRow] * (*elementThis);
    1720                }
    1721           }
    1722      }
    1723      if (fabs(value) > zeroTolerance) {
    1724           array[numberNonZero] = value;
    1725           index[numberNonZero++] = jColumn;
    1726      }
     1650  int numberThreads = abcState();
     1651  if (numberThreads) {
     1652    clpTempInfo info[ABOCA_LITE];
     1653    int chunk = (numberActiveColumns_ + numberThreads - 1) / numberThreads;
     1654    int n = 0;
     1655    for (int i = 0; i < numberThreads; i++) {
     1656      info[i].which = index + n;
     1657      info[i].infeas = array + n;
     1658      info[i].work = const_cast< double * >(pi);
     1659      info[i].status = const_cast< unsigned char * >(status);
     1660      info[i].startColumn = n;
     1661      info[i].element = elementByColumn;
     1662      info[i].start = columnStart;
     1663      info[i].row = row;
     1664      info[i].numberToDo = CoinMin(chunk, numberActiveColumns_ - n);
     1665      info[i].tolerance = zeroTolerance;
     1666      n += chunk;
     1667    }
     1668    for (int i = 0; i < numberThreads; i++) {
     1669      cilk_spawn transposeTimesUnscaledBit(info[i]);
     1670    }
     1671    cilk_sync;
     1672    for (int i = 0; i < numberThreads; i++)
     1673      numberNonZero += info[i].numberAdded;
     1674    moveAndZero(info, 2, NULL);
     1675  } else {
     1676#endif
     1677    double value = 0.0;
     1678    int jColumn = -1;
     1679    for (int iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
     1680      bool wanted = ((status[iColumn] & 3) != 1);
     1681      if (fabs(value) > zeroTolerance) {
     1682        array[numberNonZero] = value;
     1683        index[numberNonZero++] = jColumn;
     1684      }
     1685      value = 0.0;
     1686      if (wanted) {
     1687        CoinBigIndex start = columnStart[iColumn];
     1688        CoinBigIndex end = columnStart[iColumn + 1];
     1689        jColumn = iColumn;
     1690        int n = static_cast< int >(end - start);
     1691        bool odd = (n & 1) != 0;
     1692        n = n >> 1;
     1693        const int *COIN_RESTRICT rowThis = row + start;
     1694        const double *COIN_RESTRICT elementThis = elementByColumn + start;
     1695        for (; n; n--) {
     1696          int iRow0 = *rowThis;
     1697          int iRow1 = *(rowThis + 1);
     1698          rowThis += 2;
     1699          value += pi[iRow0] * (*elementThis);
     1700          value += pi[iRow1] * (*(elementThis + 1));
     1701          elementThis += 2;
     1702        }
     1703        if (odd) {
     1704          int iRow = *rowThis;
     1705          value += pi[iRow] * (*elementThis);
     1706        }
     1707      }
     1708    }
     1709    if (fabs(value) > zeroTolerance) {
     1710      array[numberNonZero] = value;
     1711      index[numberNonZero++] = jColumn;
     1712    }
    17271713#if ABOCA_LITE
    1728             }
    1729 #endif
    1730      return numberNonZero;
     1714  }
     1715#endif
     1716  return numberNonZero;
    17311717}
    17321718#if ABOCA_LITE
    17331719static void
    1734 transposeTimesUnscaledBit2(clpTempInfo & info)
    1735 {
    1736   const CoinBigIndex * COIN_RESTRICT columnStart=info.start;
    1737   const int * COIN_RESTRICT row = info.row;
    1738   const double * COIN_RESTRICT elementByColumn = info.element;
    1739   double * COIN_RESTRICT array = info.infeas;
    1740   const double * COIN_RESTRICT reducedCost = info.reducedCost;
    1741   int * COIN_RESTRICT index = info.which;
    1742   double * COIN_RESTRICT spareArray = info.spare;
    1743   int * COIN_RESTRICT spareIndex = info.index;
    1744   const unsigned char * COIN_RESTRICT status = info.status;
    1745   double zeroTolerance=info.tolerance;
    1746   double dualTolerance=info.dualTolerance;
     1720transposeTimesUnscaledBit2(clpTempInfo &info)
     1721{
     1722  const CoinBigIndex *COIN_RESTRICT columnStart = info.start;
     1723  const int *COIN_RESTRICT row = info.row;
     1724  const double *COIN_RESTRICT elementByColumn = info.element;
     1725  double *COIN_RESTRICT array = info.infeas;
     1726  const double *COIN_RESTRICT reducedCost = info.reducedCost;
     1727  int *COIN_RESTRICT index = info.which;
     1728  double *COIN_RESTRICT spareArray = info.spare;
     1729  int *COIN_RESTRICT spareIndex = info.index;
     1730  const unsigned char *COIN_RESTRICT status = info.status;
     1731  double zeroTolerance = info.tolerance;
     1732  double dualTolerance = info.dualTolerance;
    17471733  double acceptablePivot = info.acceptablePivot;
    1748   const double * COIN_RESTRICT pi = info.work;
     1734  const double *COIN_RESTRICT pi = info.work;
    17491735  double tentativeTheta = 1.0e15;
    17501736  int numberRemaining = 0;
    17511737  double upperTheta = info.upperTheta;
    17521738  int first = info.startColumn;
    1753   int last = first+info.numberToDo;
    1754   int numberNonZero=0;
    1755   double multiplier[] = { -1.0, 1.0};
    1756   double dualT = - dualTolerance;
     1739  int last = first + info.numberToDo;
     1740  int numberNonZero = 0;
     1741  double multiplier[] = { -1.0, 1.0 };
     1742  double dualT = -dualTolerance;
    17571743  for (int iColumn = first; iColumn < last; iColumn++) {
    17581744    int wanted = (status[iColumn] & 3) - 1;
     
    17601746      double value = 0.0;
    17611747      CoinBigIndex start = columnStart[iColumn];
    1762       CoinBigIndex end = columnStart[iColumn+1];
    1763       int n = static_cast<int>(end - start);
     1748      CoinBigIndex end = columnStart[iColumn + 1];
     1749      int n = static_cast< int >(end - start);
    17641750      bool odd = (n & 1) != 0;
    17651751      n = n >> 1;
    1766       const int * COIN_RESTRICT rowThis = row + start;
    1767       const double * COIN_RESTRICT elementThis = elementByColumn + start;
     1752      const int *COIN_RESTRICT rowThis = row + start;
     1753      const double *COIN_RESTRICT elementThis = elementByColumn + start;
    17681754      for (; n; n--) {
    1769         int iRow0 = *rowThis;
    1770         int iRow1 = *(rowThis + 1);
    1771         rowThis += 2;
    1772         value += pi[iRow0] * (*elementThis);
    1773         value += pi[iRow1] * (*(elementThis + 1));
    1774         elementThis += 2;
     1755        int iRow0 = *rowThis;
     1756        int iRow1 = *(rowThis + 1);
     1757        rowThis += 2;
     1758        value += pi[iRow0] * (*elementThis);
     1759        value += pi[iRow1] * (*(elementThis + 1));
     1760        elementThis += 2;
    17751761      }
    17761762      if (odd) {
    1777         int iRow = *rowThis;
    1778         value += pi[iRow] * (*elementThis);
     1763        int iRow = *rowThis;
     1764        value += pi[iRow] * (*elementThis);
    17791765      }
    17801766      if (fabs(value) > zeroTolerance) {
    1781         double mult = multiplier[wanted-1];
    1782         double alpha = value * mult;
    1783         array[numberNonZero] = value;
    1784         index[numberNonZero++] = iColumn;
    1785         if (alpha > 0.0) {
    1786           double oldValue = reducedCost[iColumn] * mult;
    1787           double value = oldValue - tentativeTheta * alpha;
    1788           if (value < dualT) {
    1789             value = oldValue - upperTheta * alpha;
    1790             if (value < dualT && alpha >= acceptablePivot) {
    1791               upperTheta = (oldValue - dualT) / alpha;
    1792               //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
    1793             }
    1794             // add to list
    1795             spareArray[numberRemaining] = alpha * mult;
    1796             spareIndex[numberRemaining++] = iColumn;
    1797           }
    1798         }
    1799       }
    1800     }
    1801   }
    1802   info.numberAdded=numberNonZero;
    1803   info.numberRemaining=numberRemaining;
    1804   info.upperTheta=upperTheta;
     1767        double mult = multiplier[wanted - 1];
     1768        double alpha = value * mult;
     1769        array[numberNonZero] = value;
     1770        index[numberNonZero++] = iColumn;
     1771        if (alpha > 0.0) {
     1772          double oldValue = reducedCost[iColumn] * mult;
     1773          double value = oldValue - tentativeTheta * alpha;
     1774          if (value < dualT) {
     1775            value = oldValue - upperTheta * alpha;
     1776            if (value < dualT && alpha >= acceptablePivot) {
     1777              upperTheta = (oldValue - dualT) / alpha;
     1778              //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
     1779            }
     1780            // add to list
     1781            spareArray[numberRemaining] = alpha * mult;
     1782            spareIndex[numberRemaining++] = iColumn;
     1783          }
     1784        }
     1785      }
     1786    }
     1787  }
     1788  info.numberAdded = numberNonZero;
     1789  info.numberRemaining = numberRemaining;
     1790  info.upperTheta = upperTheta;
    18051791}
    18061792#endif
    18071793/* Meat of transposeTimes by column when not scaled and skipping
    18081794   and doing part of dualColumn */
    1809 int
    1810 ClpPackedMatrix::gutsOfTransposeTimesUnscaled(const double * COIN_RESTRICT pi,
    1811           int * COIN_RESTRICT index,
    1812           double * COIN_RESTRICT array,
    1813           const unsigned char * COIN_RESTRICT status,
    1814           int * COIN_RESTRICT spareIndex,
    1815           double * COIN_RESTRICT spareArray,
    1816           const double * COIN_RESTRICT reducedCost,
    1817           double & upperThetaP,
    1818           double acceptablePivot,
    1819           double dualTolerance,
    1820           int & numberRemainingP,
    1821           const double zeroTolerance) const
    1822 {
    1823      int numberRemaining = numberRemainingP;
    1824      double upperTheta = upperThetaP;
    1825      int numberNonZero = 0;
    1826      // get matrix data pointers
    1827      const int * COIN_RESTRICT row = matrix_->getIndices();
    1828      const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
    1829      const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
     1795int ClpPackedMatrix::gutsOfTransposeTimesUnscaled(const double *COIN_RESTRICT pi,
     1796  int *COIN_RESTRICT index,
     1797  double *COIN_RESTRICT array,
     1798  const unsigned char *COIN_RESTRICT status,
     1799  int *COIN_RESTRICT spareIndex,
     1800  double *COIN_RESTRICT spareArray,
     1801  const double *COIN_RESTRICT reducedCost,
     1802  double &upperThetaP,
     1803  double acceptablePivot,
     1804  double dualTolerance,
     1805  int &numberRemainingP,
     1806  const double zeroTolerance) const
     1807{
     1808  int numberRemaining = numberRemainingP;
     1809  double upperTheta = upperThetaP;
     1810  int numberNonZero = 0;
     1811  // get matrix data pointers
     1812  const int *COIN_RESTRICT row = matrix_->getIndices();
     1813  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     1814  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
    18301815#if ABOCA_LITE
    1831      int numberThreads=abcState();
    1832      if (numberThreads) {
    1833      clpTempInfo info[ABOCA_LITE];
    1834      int chunk = (numberActiveColumns_+numberThreads-1)/numberThreads;
    1835      int n=0;
    1836      for (int i=0;i<numberThreads;i++) {
    1837        info[i].upperTheta = upperThetaP;
    1838        info[i].acceptablePivot = acceptablePivot;
    1839        info[i].reducedCost = const_cast<double *>(reducedCost);
    1840        info[i].index = spareIndex+n+numberRemaining;
    1841        info[i].spare = spareArray+n+numberRemaining;
    1842        info[i].which = index+n;
    1843        info[i].infeas = array+n;
    1844        info[i].work=const_cast<double *>(pi);
    1845        info[i].status=const_cast<unsigned char *>(status);
    1846        info[i].startColumn=n;
    1847        info[i].element=elementByColumn;
    1848        info[i].start=columnStart;
    1849        info[i].row=row;
    1850        info[i].numberToDo=CoinMin(chunk,numberActiveColumns_-n);
    1851        info[i].tolerance=zeroTolerance;
    1852        info[i].dualTolerance=dualTolerance;
    1853        n += chunk;
    1854      }
    1855      for (int i=0;i<numberThreads;i++) {
    1856        cilk_spawn transposeTimesUnscaledBit2(info[i]);
    1857      }
    1858      cilk_sync;
    1859      for (int i=0;i<numberThreads;i++) {
    1860        numberNonZero += info[i].numberAdded;
    1861        numberRemaining += info[i].numberRemaining;
    1862        upperTheta = CoinMin(upperTheta,static_cast<double>(info[i].upperTheta));
    1863      }
    1864      moveAndZero(info,1,NULL);
    1865      moveAndZero(info,2,NULL);
    1866      } else {
    1867 #endif
    1868      double tentativeTheta = 1.0e15;
    1869      double multiplier[] = { -1.0, 1.0};
    1870      double dualT = - dualTolerance;
    1871      for (int iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
    1872           int wanted = (status[iColumn] & 3) - 1;
    1873           if (wanted) {
    1874                double value = 0.0;
    1875                CoinBigIndex start = columnStart[iColumn];
    1876                CoinBigIndex end = columnStart[iColumn+1];
    1877                int n = static_cast<int>(end - start);
     1816  int numberThreads = abcState();
     1817  if (numberThreads) {
     1818    clpTempInfo info[ABOCA_LITE];
     1819    int chunk = (numberActiveColumns_ + numberThreads - 1) / numberThreads;
     1820    int n = 0;
     1821    for (int i = 0; i < numberThreads; i++) {
     1822      info[i].upperTheta = upperThetaP;
     1823      info[i].acceptablePivot = acceptablePivot;
     1824      info[i].reducedCost = const_cast< double * >(reducedCost);
     1825      info[i].index = spareIndex + n + numberRemaining;
     1826      info[i].spare = spareArray + n + numberRemaining;
     1827      info[i].which = index + n;
     1828      info[i].infeas = array + n;
     1829      info[i].work = const_cast< double * >(pi);
     1830      info[i].status = const_cast< unsigned char * >(status);
     1831      info[i].startColumn = n;
     1832      info[i].element = elementByColumn;
     1833      info[i].start = columnStart;
     1834      info[i].row = row;
     1835      info[i].numberToDo = CoinMin(chunk, numberActiveColumns_ - n);
     1836      info[i].tolerance = zeroTolerance;
     1837      info[i].dualTolerance = dualTolerance;
     1838      n += chunk;
     1839    }
     1840    for (int i = 0; i < numberThreads; i++) {
     1841      cilk_spawn transposeTimesUnscaledBit2(info[i]);
     1842    }
     1843    cilk_sync;
     1844    for (int i = 0; i < numberThreads; i++) {
     1845      numberNonZero += info[i].numberAdded;
     1846      numberRemaining += info[i].numberRemaining;
     1847      upperTheta = CoinMin(upperTheta, static_cast< double >(info[i].upperTheta));
     1848    }
     1849    moveAndZero(info, 1, NULL);
     1850    moveAndZero(info, 2, NULL);
     1851  } else {
     1852#endif
     1853    double tentativeTheta = 1.0e15;
     1854    double multiplier[] = { -1.0, 1.0 };
     1855    double dualT = -dualTolerance;
     1856    for (int iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
     1857      int wanted = (status[iColumn] & 3) - 1;
     1858      if (wanted) {
     1859        double value = 0.0;
     1860        CoinBigIndex start = columnStart[iColumn];
     1861        CoinBigIndex end = columnStart[iColumn + 1];
     1862        int n = static_cast< int >(end - start);
    18781863#if 1
    1879                bool odd = (n & 1) != 0;
    1880                n = n >> 1;
    1881                const int * COIN_RESTRICT rowThis = row + start;
    1882                const double * COIN_RESTRICT elementThis = elementByColumn + start;
    1883                for (; n; n--) {
    1884                     int iRow0 = *rowThis;
    1885                     int iRow1 = *(rowThis + 1);
    1886                     rowThis += 2;
    1887                     value += pi[iRow0] * (*elementThis);
    1888                     value += pi[iRow1] * (*(elementThis + 1));
    1889                     elementThis += 2;
    1890                }
    1891                if (odd) {
    1892                     int iRow = *rowThis;
    1893                     value += pi[iRow] * (*elementThis);
    1894                }
     1864        bool odd = (n & 1) != 0;
     1865        n = n >> 1;
     1866        const int *COIN_RESTRICT rowThis = row + start;
     1867        const double *COIN_RESTRICT elementThis = elementByColumn + start;
     1868        for (; n; n--) {
     1869          int iRow0 = *rowThis;
     1870          int iRow1 = *(rowThis + 1);
     1871          rowThis += 2;
     1872          value += pi[iRow0] * (*elementThis);
     1873          value += pi[iRow1] * (*(elementThis + 1));
     1874          elementThis += 2;
     1875        }
     1876        if (odd) {
     1877          int iRow = *rowThis;
     1878          value += pi[iRow] * (*elementThis);
     1879        }
    18951880#else
    1896                const int * COIN_RESTRICT rowThis = &row[end-16];
    1897                const double * COIN_RESTRICT elementThis = &elementByColumn[end-16];
    1898                bool odd = (n & 1) != 0;
    1899                n = n >> 1;
    1900                double value2 = 0.0;
    1901                if (odd) {
    1902                     int iRow = row[start];
    1903                     value2 = pi[iRow] * elementByColumn[start];
    1904                }
    1905                switch (n) {
    1906                default: {
    1907                     if (odd) {
    1908                          start++;
    1909                     }
    1910                     n -= 8;
    1911                     for (; n; n--) {
    1912                          int iRow0 = row[start];
    1913                          int iRow1 = row[start+1];
    1914                          value += pi[iRow0] * elementByColumn[start];
    1915                          value2 += pi[iRow1] * elementByColumn[start+1];
    1916                          start += 2;
    1917                     }
    1918                     case 8: {
    1919                          int iRow0 = rowThis[16-16];
    1920                          int iRow1 = rowThis[16-15];
    1921                          value += pi[iRow0] * elementThis[16-16];
    1922                          value2 += pi[iRow1] * elementThis[16-15];
    1923                     }
    1924                     case 7: {
    1925                          int iRow0 = rowThis[16-14];
    1926                          int iRow1 = rowThis[16-13];
    1927                          value += pi[iRow0] * elementThis[16-14];
    1928                          value2 += pi[iRow1] * elementThis[16-13];
    1929                     }
    1930                     case 6: {
    1931                          int iRow0 = rowThis[16-12];
    1932                          int iRow1 = rowThis[16-11];
    1933                          value += pi[iRow0] * elementThis[16-12];
    1934                          value2 += pi[iRow1] * elementThis[16-11];
    1935                     }
    1936                     case 5: {
    1937                          int iRow0 = rowThis[16-10];
    1938                          int iRow1 = rowThis[16-9];
    1939                          value += pi[iRow0] * elementThis[16-10];
    1940                          value2 += pi[iRow1] * elementThis[16-9];
    1941                     }
    1942                     case 4: {
    1943                          int iRow0 = rowThis[16-8];
    1944                          int iRow1 = rowThis[16-7];
    1945                          value += pi[iRow0] * elementThis[16-8];
    1946                          value2 += pi[iRow1] * elementThis[16-7];
    1947                     }
    1948                     case 3: {
    1949                          int iRow0 = rowThis[16-6];
    1950                          int iRow1 = rowThis[16-5];
    1951                          value += pi[iRow0] * elementThis[16-6];
    1952                          value2 += pi[iRow1] * elementThis[16-5];
    1953                     }
    1954                     case 2: {
    1955                          int iRow0 = rowThis[16-4];
    1956                          int iRow1 = rowThis[16-3];
    1957                          value += pi[iRow0] * elementThis[16-4];
    1958                          value2 += pi[iRow1] * elementThis[16-3];
    1959                     }
    1960                     case 1: {
    1961                          int iRow0 = rowThis[16-2];
    1962                          int iRow1 = rowThis[16-1];
    1963                          value += pi[iRow0] * elementThis[16-2];
    1964                          value2 += pi[iRow1] * elementThis[16-1];
    1965                     }
    1966                     case 0:
    1967                          ;
    1968                     }
    1969                }
    1970                value += value2;
    1971 #endif
    1972                if (fabs(value) > zeroTolerance) {
    1973                     double mult = multiplier[wanted-1];
    1974                     double alpha = value * mult;
    1975                     array[numberNonZero] = value;
    1976                     index[numberNonZero++] = iColumn;
    1977                     if (alpha > 0.0) {
    1978                          double oldValue = reducedCost[iColumn] * mult;
    1979                          double value = oldValue - tentativeTheta * alpha;
    1980                          if (value < dualT) {
    1981                               value = oldValue - upperTheta * alpha;
    1982                               if (value < dualT && alpha >= acceptablePivot) {
    1983                                    upperTheta = (oldValue - dualT) / alpha;
    1984                                    //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
    1985                               }
    1986                               // add to list
    1987                               spareArray[numberRemaining] = alpha * mult;
    1988                               spareIndex[numberRemaining++] = iColumn;
    1989                          }
    1990                     }
    1991                }
    1992           }
    1993      }
     1881      const int *COIN_RESTRICT rowThis = &row[end - 16];
     1882      const double *COIN_RESTRICT elementThis = &elementByColumn[end - 16];
     1883      bool odd = (n & 1) != 0;
     1884      n = n >> 1;
     1885      double value2 = 0.0;
     1886      if (odd) {
     1887        int iRow = row[start];
     1888        value2 = pi[iRow] * elementByColumn[start];
     1889      }
     1890      switch (n) {
     1891      default: {
     1892        if (odd) {
     1893          start++;
     1894        }
     1895        n -= 8;
     1896        for (; n; n--) {
     1897          int iRow0 = row[start];
     1898          int iRow1 = row[start + 1];
     1899          value += pi[iRow0] * elementByColumn[start];
     1900          value2 += pi[iRow1] * elementByColumn[start + 1];
     1901          start += 2;
     1902        }
     1903      case 8: {
     1904        int iRow0 = rowThis[16 - 16];
     1905        int iRow1 = rowThis[16 - 15];
     1906        value += pi[iRow0] * elementThis[16 - 16];
     1907        value2 += pi[iRow1] * elementThis[16 - 15];
     1908      }
     1909      case 7: {
     1910        int iRow0 = rowThis[16 - 14];
     1911        int iRow1 = rowThis[16 - 13];
     1912        value += pi[iRow0] * elementThis[16 - 14];
     1913        value2 += pi[iRow1] * elementThis[16 - 13];
     1914      }
     1915      case 6: {
     1916        int iRow0 = rowThis[16 - 12];
     1917        int iRow1 = rowThis[16 - 11];
     1918        value += pi[iRow0] * elementThis[16 - 12];
     1919        value2 += pi[iRow1] * elementThis[16 - 11];
     1920      }
     1921      case 5: {
     1922        int iRow0 = rowThis[16 - 10];
     1923        int iRow1 = rowThis[16 - 9];
     1924        value += pi[iRow0] * elementThis[16 - 10];
     1925        value2 += pi[iRow1] * elementThis[16 - 9];
     1926      }
     1927      case 4: {
     1928        int iRow0 = rowThis[16 - 8];
     1929        int iRow1 = rowThis[16 - 7];
     1930        value += pi[iRow0] * elementThis[16 - 8];
     1931        value2 += pi[iRow1] * elementThis[16 - 7];
     1932      }
     1933      case 3: {
     1934        int iRow0 = rowThis[16 - 6];
     1935        int iRow1 = rowThis[16 - 5];
     1936        value += pi[iRow0] * elementThis[16 - 6];
     1937        value2 += pi[iRow1] * elementThis[16 - 5];
     1938      }
     1939      case 2: {
     1940        int iRow0 = rowThis[16 - 4];
     1941        int iRow1 = rowThis[16 - 3];
     1942        value += pi[iRow0] * elementThis[16 - 4];
     1943        value2 += pi[iRow1] * elementThis[16 - 3];
     1944      }
     1945      case 1: {
     1946        int iRow0 = rowThis[16 - 2];
     1947        int iRow1 = rowThis[16 - 1];
     1948        value += pi[iRow0] * elementThis[16 - 2];
     1949        value2 += pi[iRow1] * elementThis[16 - 1];
     1950      }
     1951      case 0:;
     1952      }
     1953      }
     1954      value += value2;
     1955#endif
     1956        if (fabs(value) > zeroTolerance) {
     1957          double mult = multiplier[wanted - 1];
     1958          double alpha = value * mult;
     1959          array[numberNonZero] = value;
     1960          index[numberNonZero++] = iColumn;
     1961          if (alpha > 0.0) {
     1962            double oldValue = reducedCost[iColumn] * mult;
     1963            double value = oldValue - tentativeTheta * alpha;
     1964            if (value < dualT) {
     1965              value = oldValue - upperTheta * alpha;
     1966              if (value < dualT && alpha >= acceptablePivot) {
     1967                upperTheta = (oldValue - dualT) / alpha;
     1968                //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
     1969              }
     1970              // add to list
     1971              spareArray[numberRemaining] = alpha * mult;
     1972              spareIndex[numberRemaining++] = iColumn;
     1973            }
     1974          }
     1975        }
     1976      }
     1977    }
    19941978#if ABOCA_LITE
    1995             }
    1996 #endif
    1997      numberRemainingP = numberRemaining;
    1998      upperThetaP = upperTheta;
    1999      return numberNonZero;
     1979  }
     1980#endif
     1981  numberRemainingP = numberRemaining;
     1982  upperThetaP = upperTheta;
     1983  return numberNonZero;
    20001984}
    20011985// Meat of transposeTimes by column when scaled
    2002 int
    2003 ClpPackedMatrix::gutsOfTransposeTimesScaled(const double * COIN_RESTRICT pi,
    2004           const double * COIN_RESTRICT columnScale,
    2005           int * COIN_RESTRICT index,
    2006           double * COIN_RESTRICT array,
    2007           const unsigned char * COIN_RESTRICT status,                            const double zeroTolerance) const
    2008 {
    2009      int numberNonZero = 0;
    2010      // get matrix data pointers
    2011      const int * COIN_RESTRICT row = matrix_->getIndices();
    2012      const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
    2013      const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
    2014      double value = 0.0;
    2015      int jColumn = -1;
    2016      for (int iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
    2017           bool wanted = ((status[iColumn] & 3) != 1);
    2018           if (fabs(value) > zeroTolerance) {
    2019                array[numberNonZero] = value;
    2020                index[numberNonZero++] = jColumn;
    2021           }
    2022           value = 0.0;
    2023           if (wanted) {
    2024                double scale = columnScale[iColumn];
    2025                CoinBigIndex start = columnStart[iColumn];
    2026                CoinBigIndex end = columnStart[iColumn+1];
    2027                jColumn = iColumn;
    2028                for (CoinBigIndex j = start; j < end; j++) {
    2029                     int iRow = row[j];
    2030                     value += pi[iRow] * elementByColumn[j];
    2031                }
    2032                value *= scale;
    2033           }
    2034      }
    2035      if (fabs(value) > zeroTolerance) {
    2036           array[numberNonZero] = value;
    2037           index[numberNonZero++] = jColumn;
    2038      }
    2039      return numberNonZero;
     1986int ClpPackedMatrix::gutsOfTransposeTimesScaled(const double *COIN_RESTRICT pi,
     1987  const double *COIN_RESTRICT columnScale,
     1988  int *COIN_RESTRICT index,
     1989  double *COIN_RESTRICT array,
     1990  const unsigned char *COIN_RESTRICT status, const double zeroTolerance) const
     1991{
     1992  int numberNonZero = 0;
     1993  // get matrix data pointers
     1994  const int *COIN_RESTRICT row = matrix_->getIndices();
     1995  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     1996  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
     1997  double value = 0.0;
     1998  int jColumn = -1;
     1999  for (int iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
     2000    bool wanted = ((status[iColumn] & 3) != 1);
     2001    if (fabs(value) > zeroTolerance) {
     2002      array[numberNonZero] = value;
     2003      index[numberNonZero++] = jColumn;
     2004    }
     2005    value = 0.0;
     2006    if (wanted) {
     2007      double scale = columnScale[iColumn];
     2008      CoinBigIndex start = columnStart[iColumn];
     2009      CoinBigIndex end = columnStart[iColumn + 1];
     2010      jColumn = iColumn;
     2011      for (CoinBigIndex j = start; j < end; j++) {
     2012        int iRow = row[j];
     2013        value += pi[iRow] * elementByColumn[j];
     2014      }
     2015      value *= scale;
     2016    }
     2017  }
     2018  if (fabs(value) > zeroTolerance) {
     2019    array[numberNonZero] = value;
     2020    index[numberNonZero++] = jColumn;
     2021  }
     2022  return numberNonZero;
    20402023}
    20412024#if ABOCA_LITE
    20422025static void
    2043 packDownBit(clpTempInfo & info)
    2044 {
    2045   double * COIN_RESTRICT output = info.infeas;
    2046   int * COIN_RESTRICT index = info.which;
    2047   double zeroTolerance=info.tolerance;
     2026packDownBit(clpTempInfo &info)
     2027{
     2028  double *COIN_RESTRICT output = info.infeas;
     2029  int *COIN_RESTRICT index = info.which;
     2030  double zeroTolerance = info.tolerance;
    20482031  int first = info.startColumn;
    20492032  int last = info.numberToDo;
    2050   int numberNonZero=0;
     2033  int numberNonZero = 0;
    20512034  for (int i = 0; i < last; i++) {
    20522035    double value = output[i];
     
    20542037      output[i] = 0.0;
    20552038      if (fabs(value) > zeroTolerance) {
    2056         output[numberNonZero] = value;
    2057         index[numberNonZero++] = i+first;
    2058       }
    2059     }
    2060   }
    2061   info.numberAdded=numberNonZero;
     2039        output[numberNonZero] = value;
     2040        index[numberNonZero++] = i + first;
     2041      }
     2042    }
     2043  }
     2044  info.numberAdded = numberNonZero;
    20622045}
    20632046#endif
    20642047// Meat of transposeTimes by row n > K if packed - returns number nonzero
    2065 int
    2066 ClpPackedMatrix::gutsOfTransposeTimesByRowGEK(const CoinIndexedVector * COIN_RESTRICT piVector,
    2067           int * COIN_RESTRICT index,
    2068           double * COIN_RESTRICT output,
    2069           int numberColumns,
    2070           const double tolerance,
    2071           const double scalar) const
    2072 {
    2073      const double * COIN_RESTRICT pi = piVector->denseVector();
    2074      int numberInRowArray = piVector->getNumElements();
    2075      const int * COIN_RESTRICT column = matrix_->getIndices();
    2076      const CoinBigIndex * COIN_RESTRICT rowStart = matrix_->getVectorStarts();
    2077      const double * COIN_RESTRICT element = matrix_->getElements();
    2078      const int * COIN_RESTRICT whichRow = piVector->getIndices();
    2079      // ** Row copy is already scaled
    2080      for (int i = 0; i < numberInRowArray; i++) {
    2081           int iRow = whichRow[i];
    2082           double value = pi[i] * scalar;
    2083           CoinBigIndex start = rowStart[iRow];
    2084           CoinBigIndex end = rowStart[iRow+1];
    2085           int n = static_cast<int>(end - start);
    2086           const int * COIN_RESTRICT columnThis = column + start;
    2087           const double * COIN_RESTRICT elementThis = element + start;
     2048int ClpPackedMatrix::gutsOfTransposeTimesByRowGEK(const CoinIndexedVector *COIN_RESTRICT piVector,
     2049  int *COIN_RESTRICT index,
     2050  double *COIN_RESTRICT output,
     2051  int numberColumns,
     2052  const double tolerance,
     2053  const double scalar) const
     2054{
     2055  const double *COIN_RESTRICT pi = piVector->denseVector();
     2056  int numberInRowArray = piVector->getNumElements();
     2057  const int *COIN_RESTRICT column = matrix_->getIndices();
     2058  const CoinBigIndex *COIN_RESTRICT rowStart = matrix_->getVectorStarts();
     2059  const double *COIN_RESTRICT element = matrix_->getElements();
     2060  const int *COIN_RESTRICT whichRow = piVector->getIndices();
     2061  // ** Row copy is already scaled
     2062  for (int i = 0; i < numberInRowArray; i++) {
     2063    int iRow = whichRow[i];
     2064    double value = pi[i] * scalar;
     2065    CoinBigIndex start = rowStart[iRow];
     2066    CoinBigIndex end = rowStart[iRow + 1];
     2067    int n = static_cast< int >(end - start);
     2068    const int *COIN_RESTRICT columnThis = column + start;
     2069    const double *COIN_RESTRICT elementThis = element + start;
    20882070
    2089           // could do by twos
    2090           for (; n; n--) {
    2091                int iColumn = *columnThis;
    2092                columnThis++;
    2093                double elValue = *elementThis;
    2094                elementThis++;
    2095                elValue *= value;
    2096                output[iColumn] += elValue;
    2097           }
    2098      }
    2099      // get rid of tiny values and count
    2100      int numberNonZero = 0;
     2071    // could do by twos
     2072    for (; n; n--) {
     2073      int iColumn = *columnThis;
     2074      columnThis++;
     2075      double elValue = *elementThis;
     2076      elementThis++;
     2077      elValue *= value;
     2078      output[iColumn] += elValue;
     2079    }
     2080  }
     2081  // get rid of tiny values and count
     2082  int numberNonZero = 0;
    21012083#if ABOCA_LITE
    2102      int numberThreads=abcState();
    2103      if (numberThreads) {
    2104      clpTempInfo info[ABOCA_LITE];
    2105      int chunk = (numberColumns+numberThreads-1)/numberThreads;
    2106      int n=0;
    2107      for (int i=0;i<numberThreads;i++) {
    2108        info[i].which = index+n;
    2109        info[i].infeas = output+n;
    2110        info[i].startColumn=n;
    2111        info[i].numberToDo=CoinMin(chunk,numberColumns-n);
    2112        info[i].tolerance=tolerance;
    2113        n += chunk;
    2114      }
    2115      for (int i=0;i<numberThreads;i++) {
    2116        cilk_spawn packDownBit(info[i]);
    2117      }
    2118      cilk_sync;
    2119      for (int i=0;i<numberThreads;i++) {
    2120        numberNonZero += info[i].numberAdded;
    2121      }
    2122      moveAndZero(info,2,NULL);
    2123      } else {
    2124 #endif
    2125      for (int i = 0; i < numberColumns; i++) {
    2126           double value = output[i];
    2127           if (value) {
    2128                output[i] = 0.0;
    2129                if (fabs(value) > tolerance) {
    2130                     output[numberNonZero] = value;
    2131                     index[numberNonZero++] = i;
    2132                }
    2133           }
    2134      }
     2084  int numberThreads = abcState();
     2085  if (numberThreads) {
     2086    clpTempInfo info[ABOCA_LITE];
     2087    int chunk = (numberColumns + numberThreads - 1) / numberThreads;
     2088    int n = 0;
     2089    for (int i = 0; i < numberThreads; i++) {
     2090      info[i].which = index + n;
     2091      info[i].infeas = output + n;
     2092      info[i].startColumn = n;
     2093      info[i].numberToDo = CoinMin(chunk, numberColumns - n);
     2094      info[i].tolerance = tolerance;
     2095      n += chunk;
     2096    }
     2097    for (int i = 0; i < numberThreads; i++) {
     2098      cilk_spawn packDownBit(info[i]);
     2099    }
     2100    cilk_sync;
     2101    for (int i = 0; i < numberThreads; i++) {
     2102      numberNonZero += info[i].numberAdded;
     2103    }
     2104    moveAndZero(info, 2, NULL);
     2105  } else {
     2106#endif
     2107    for (int i = 0; i < numberColumns; i++) {
     2108      double value = output[i];
     2109      if (value) {
     2110        output[i] = 0.0;
     2111        if (fabs(value) > tolerance) {
     2112          output[numberNonZero] = value;
     2113          index[numberNonZero++] = i;
     2114        }
     2115      }
     2116    }
    21352117#if ABOCA_LITE
    2136             }
     2118  }
    21372119#endif
    21382120#ifndef NDEBUG
    2139      for (int i = numberNonZero; i < numberColumns; i++)
    2140           assert(!output[i]);
    2141 #endif
    2142      return numberNonZero;
     2121  for (int i = numberNonZero; i < numberColumns; i++)
     2122    assert(!output[i]);
     2123#endif
     2124  return numberNonZero;
    21432125}
    21442126// Meat of transposeTimes by row n == 2 if packed
    2145 void
    2146 ClpPackedMatrix::gutsOfTransposeTimesByRowEQ2(const CoinIndexedVector * piVector, CoinIndexedVector * output,
    2147           CoinIndexedVector * spareVector, const double tolerance, const double scalar) const
    2148 {
    2149      double * COIN_RESTRICT pi = piVector->denseVector();
    2150      int numberNonZero = 0;
    2151      int * COIN_RESTRICT index = output->getIndices();
    2152      double * COIN_RESTRICT array = output->denseVector();
    2153      const int * COIN_RESTRICT column = matrix_->getIndices();
    2154      const CoinBigIndex * COIN_RESTRICT rowStart = matrix_->getVectorStarts();
    2155      const double * COIN_RESTRICT element = matrix_->getElements();
    2156      const int * COIN_RESTRICT whichRow = piVector->getIndices();
    2157      int iRow0 = whichRow[0];
    2158      int iRow1 = whichRow[1];
    2159      double pi0 = pi[0];
    2160      double pi1 = pi[1];
    2161      if (rowStart[iRow0+1] - rowStart[iRow0] >
    2162                rowStart[iRow1+1] - rowStart[iRow1]) {
    2163           // do one with fewer first
    2164           iRow0 = iRow1;
    2165           iRow1 = whichRow[0];
    2166           pi0 = pi1;
    2167           pi1 = pi[0];
    2168      }
    2169      // and set up mark as char array
    2170      char * COIN_RESTRICT marked = reinterpret_cast<char *> (index + output->capacity());
    2171      int * COIN_RESTRICT lookup = spareVector->getIndices();
    2172      double value = pi0 * scalar;
    2173      CoinBigIndex j;
    2174      for (j = rowStart[iRow0]; j < rowStart[iRow0+1]; j++) {
    2175           int iColumn = column[j];
    2176           double elValue = element[j];
    2177           double value2 = value * elValue;
    2178           array[numberNonZero] = value2;
    2179           marked[iColumn] = 1;
    2180           lookup[iColumn] = numberNonZero;
    2181           index[numberNonZero++] = iColumn;
    2182      }
    2183      value = pi1 * scalar;
    2184      for (j = rowStart[iRow1]; j < rowStart[iRow1+1]; j++) {
    2185           int iColumn = column[j];
    2186           double elValue = element[j];
    2187           double value2 = value * elValue;
    2188           // I am assuming no zeros in matrix
    2189           if (marked[iColumn]) {
    2190                int iLookup = lookup[iColumn];
    2191                array[iLookup] += value2;
    2192           } else {
    2193                if (fabs(value2) > tolerance) {
    2194                     array[numberNonZero] = value2;
    2195                     index[numberNonZero++] = iColumn;
    2196                }
    2197           }
    2198      }
    2199      // get rid of tiny values and zero out marked
    2200      int i;
    2201      int n = numberNonZero;
    2202      numberNonZero=0;
    2203      for (i = 0; i < n; i++) {
    2204           int iColumn = index[i];
    2205           marked[iColumn] = 0;
    2206           if (fabs(array[i]) > tolerance) {
    2207             array[numberNonZero] = array[i];
    2208             index[numberNonZero++] = index[i];
    2209           }
    2210      }
    2211      memset(array+numberNonZero,0,(n-numberNonZero)*sizeof(double));
    2212      output->setNumElements(numberNonZero);
    2213      spareVector->setNumElements(0);
     2127void ClpPackedMatrix::gutsOfTransposeTimesByRowEQ2(const CoinIndexedVector *piVector, CoinIndexedVector *output,
     2128  CoinIndexedVector *spareVector, const double tolerance, const double scalar) const
     2129{
     2130  double *COIN_RESTRICT pi = piVector->denseVector();
     2131  int numberNonZero = 0;
     2132  int *COIN_RESTRICT index = output->getIndices();
     2133  double *COIN_RESTRICT array = output->denseVector();
     2134  const int *COIN_RESTRICT column = matrix_->getIndices();
     2135  const CoinBigIndex *COIN_RESTRICT rowStart = matrix_->getVectorStarts();
     2136  const double *COIN_RESTRICT element = matrix_->getElements();
     2137  const int *COIN_RESTRICT whichRow = piVector->getIndices();
     2138  int iRow0 = whichRow[0];
     2139  int iRow1 = whichRow[1];
     2140  double pi0 = pi[0];
     2141  double pi1 = pi[1];
     2142  if (rowStart[iRow0 + 1] - rowStart[iRow0] > rowStart[iRow1 + 1] - rowStart[iRow1]) {
     2143    // do one with fewer first
     2144    iRow0 = iRow1;
     2145    iRow1 = whichRow[0];
     2146    pi0 = pi1;
     2147    pi1 = pi[0];
     2148  }
     2149  // and set up mark as char array
     2150  char *COIN_RESTRICT marked = reinterpret_cast< char * >(index + output->capacity());
     2151  int *COIN_RESTRICT lookup = spareVector->getIndices();
     2152  double value = pi0 * scalar;
     2153  CoinBigIndex j;
     2154  for (j = rowStart[iRow0]; j < rowStart[iRow0 + 1]; j++) {
     2155    int iColumn = column[j];
     2156    double elValue = element[j];
     2157    double value2 = value * elValue;
     2158    array[numberNonZero] = value2;
     2159    marked[iColumn] = 1;
     2160    lookup[iColumn] = numberNonZero;
     2161    index[numberNonZero++] = iColumn;
     2162  }
     2163  value = pi1 * scalar;
     2164  for (j = rowStart[iRow1]; j < rowStart[iRow1 + 1]; j++) {
     2165    int iColumn = column[j];
     2166    double elValue = element[j];
     2167    double value2 = value * elValue;
     2168    // I am assuming no zeros in matrix
     2169    if (marked[iColumn]) {
     2170      int iLookup = lookup[iColumn];
     2171      array[iLookup] += value2;
     2172    } else {
     2173      if (fabs(value2) > tolerance) {
     2174        array[numberNonZero] = value2;
     2175        index[numberNonZero++] = iColumn;
     2176      }
     2177    }
     2178  }
     2179  // get rid of tiny values and zero out marked
     2180  int i;
     2181  int n = numberNonZero;
     2182  numberNonZero = 0;
     2183  for (i = 0; i < n; i++) {
     2184    int iColumn = index[i];
     2185    marked[iColumn] = 0;
     2186    if (fabs(array[i]) > tolerance) {
     2187      array[numberNonZero] = array[i];
     2188      index[numberNonZero++] = index[i];
     2189    }
     2190  }
     2191  memset(array + numberNonZero, 0, (n - numberNonZero) * sizeof(double));
     2192  output->setNumElements(numberNonZero);
     2193  spareVector->setNumElements(0);
    22142194}
    22152195// Meat of transposeTimes by row n == 1 if packed
    2216 void
    2217 ClpPackedMatrix::gutsOfTransposeTimesByRowEQ1(const CoinIndexedVector * piVector, CoinIndexedVector * output,
    2218           const double tolerance, const double scalar) const
    2219 {
    2220      double * COIN_RESTRICT pi = piVector->denseVector();
    2221      int numberNonZero = 0;
    2222      int * COIN_RESTRICT index = output->getIndices();
    2223      double * COIN_RESTRICT array = output->denseVector();
    2224      const int * COIN_RESTRICT column = matrix_->getIndices();
    2225      const CoinBigIndex * COIN_RESTRICT rowStart = matrix_->getVectorStarts();
    2226      const double * COIN_RESTRICT element = matrix_->getElements();
    2227      int iRow = piVector->getIndices()[0];
    2228      numberNonZero = 0;
    2229      CoinBigIndex j;
    2230      double value = pi[0] * scalar;
    2231      for (j = rowStart[iRow]; j < rowStart[iRow+1]; j++) {
    2232           int iColumn = column[j];
    2233           double elValue = element[j];
    2234           double value2 = value * elValue;
    2235           if (fabs(value2) > tolerance) {
    2236                array[numberNonZero] = value2;
    2237                index[numberNonZero++] = iColumn;
    2238           }
    2239      }
    2240      output->setNumElements(numberNonZero);
     2196void ClpPackedMatrix::gutsOfTransposeTimesByRowEQ1(const CoinIndexedVector *piVector, CoinIndexedVector *output,
     2197  const double tolerance, const double scalar) const
     2198{
     2199  double *COIN_RESTRICT pi = piVector->denseVector();
     2200  int numberNonZero = 0;
     2201  int *COIN_RESTRICT index = output->getIndices();
     2202  double *COIN_RESTRICT array = output->denseVector();
     2203  const int *COIN_RESTRICT column = matrix_->getIndices();
     2204  const CoinBigIndex *COIN_RESTRICT rowStart = matrix_->getVectorStarts();
     2205  const double *COIN_RESTRICT element = matrix_->getElements();
     2206  int iRow = piVector->getIndices()[0];
     2207  numberNonZero = 0;
     2208  CoinBigIndex j;
     2209  double value = pi[0] * scalar;
     2210  for (j = rowStart[iRow]; j < rowStart[iRow + 1]; j++) {
     2211    int iColumn = column[j];
     2212    double elValue = element[j];
     2213    double value2 = value * elValue;
     2214    if (fabs(value2) > tolerance) {
     2215      array[numberNonZero] = value2;
     2216      index[numberNonZero++] = iColumn;
     2217    }
     2218  }
     2219  output->setNumElements(numberNonZero);
    22412220}
    22422221/* Return <code>x *A in <code>z</code> but
    22432222   just for indices in y.
    22442223   Squashes small elements and knows about ClpSimplex */
    2245 void
    2246 ClpPackedMatrix::subsetTransposeTimes(const ClpSimplex * model,
    2247                                       const CoinIndexedVector * rowArray,
    2248                                       const CoinIndexedVector * y,
    2249                                       CoinIndexedVector * columnArray) const
    2250 {
    2251      columnArray->clear();
    2252      double * COIN_RESTRICT pi = rowArray->denseVector();
    2253      double * COIN_RESTRICT array = columnArray->denseVector();
    2254      int jColumn;
    2255      // get matrix data pointers
    2256      const int * COIN_RESTRICT row = matrix_->getIndices();
    2257      const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
    2258      const int * COIN_RESTRICT columnLength = matrix_->getVectorLengths();
    2259      const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
    2260      const double * COIN_RESTRICT rowScale = model->rowScale();
    2261      int numberToDo = y->getNumElements();
    2262      const int * COIN_RESTRICT which = y->getIndices();
    2263      assert (!rowArray->packedMode());
    2264      columnArray->setPacked();
    2265      ClpPackedMatrix * scaledMatrix = model->clpScaledMatrix();
    2266      int flags = flags_;
    2267      if (rowScale && scaledMatrix && !(scaledMatrix->flags() & 2)) {
    2268           flags = 0;
    2269           rowScale = NULL;
    2270           // get matrix data pointers
    2271           row = scaledMatrix->getIndices();
    2272           columnStart = scaledMatrix->getVectorStarts();
    2273           elementByColumn = scaledMatrix->getElements();
    2274      }
    2275      if (!(flags & 2) && numberToDo>2) {
    2276           // no gaps
    2277           if (!rowScale) {
    2278                int iColumn = which[0];
    2279                double value = 0.0;
    2280                CoinBigIndex j;
    2281                int columnNext = which[1];
    2282                CoinBigIndex startNext=columnStart[columnNext];
    2283                //coin_prefetch_const(row+startNext);
    2284                //coin_prefetch_const(elementByColumn+startNext);
    2285                CoinBigIndex endNext=columnStart[columnNext+1];
    2286                for (j = columnStart[iColumn];
    2287                          j < columnStart[iColumn+1]; j++) {
    2288                     int iRow = row[j];
    2289                     value += pi[iRow] * elementByColumn[j];
    2290                }
    2291                for (jColumn = 0; jColumn < numberToDo - 2; jColumn++) {
    2292                     CoinBigIndex start = startNext;
    2293                     CoinBigIndex end = endNext;
    2294                     columnNext = which[jColumn+2];
    2295                     startNext=columnStart[columnNext];
    2296                     //coin_prefetch_const(row+startNext);
    2297                     //coin_prefetch_const(elementByColumn+startNext);
    2298                     endNext=columnStart[columnNext+1];
    2299                     array[jColumn] = value;
    2300                     value = 0.0;
    2301                     for (j = start; j < end; j++) {
    2302                          int iRow = row[j];
    2303                          value += pi[iRow] * elementByColumn[j];
    2304                     }
    2305                }
    2306                array[jColumn++] = value;
    2307                value = 0.0;
    2308                for (j = startNext; j < endNext; j++) {
    2309                  int iRow = row[j];
    2310                  value += pi[iRow] * elementByColumn[j];
    2311                }
    2312                array[jColumn] = value;
    2313           } else {
     2224void ClpPackedMatrix::subsetTransposeTimes(const ClpSimplex *model,
     2225  const CoinIndexedVector *rowArray,
     2226  const CoinIndexedVector *y,
     2227  CoinIndexedVector *columnArray) const
     2228{
     2229  columnArray->clear();
     2230  double *COIN_RESTRICT pi = rowArray->denseVector();
     2231  double *COIN_RESTRICT array = columnArray->denseVector();
     2232  int jColumn;
     2233  // get matrix data pointers
     2234  const int *COIN_RESTRICT row = matrix_->getIndices();
     2235  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     2236  const int *COIN_RESTRICT columnLength = matrix_->getVectorLengths();
     2237  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
     2238  const double *COIN_RESTRICT rowScale = model->rowScale();
     2239  int numberToDo = y->getNumElements();
     2240  const int *COIN_RESTRICT which = y->getIndices();
     2241  assert(!rowArray->packedMode());
     2242  columnArray->setPacked();
     2243  ClpPackedMatrix *scaledMatrix = model->clpScaledMatrix();
     2244  int flags = flags_;
     2245  if (rowScale && scaledMatrix && !(scaledMatrix->flags() & 2)) {
     2246    flags = 0;
     2247    rowScale = NULL;
     2248    // get matrix data pointers
     2249    row = scaledMatrix->getIndices();
     2250    columnStart = scaledMatrix->getVectorStarts();
     2251    elementByColumn = scaledMatrix->getElements();
     2252  }
     2253  if (!(flags & 2) && numberToDo > 2) {
     2254    // no gaps
     2255    if (!rowScale) {
     2256      int iColumn = which[0];
     2257      double value = 0.0;
     2258      CoinBigIndex j;
     2259      int columnNext = which[1];
     2260      CoinBigIndex startNext = columnStart[columnNext];
     2261      //coin_prefetch_const(row+startNext);
     2262      //coin_prefetch_const(elementByColumn+startNext);
     2263      CoinBigIndex endNext = columnStart[columnNext + 1];
     2264      for (j = columnStart[iColumn];
     2265           j < columnStart[iColumn + 1]; j++) {
     2266        int iRow = row[j];
     2267        value += pi[iRow] * elementByColumn[j];
     2268      }
     2269      for (jColumn = 0; jColumn < numberToDo - 2; jColumn++) {
     2270        CoinBigIndex start = startNext;
     2271        CoinBigIndex end = endNext;
     2272        columnNext = which[jColumn + 2];
     2273        startNext = columnStart[columnNext];
     2274        //coin_prefetch_const(row+startNext);
     2275        //coin_prefetch_const(elementByColumn+startNext);
     2276        endNext = columnStart[columnNext + 1];
     2277        array[jColumn] = value;
     2278        value = 0.0;
     2279        for (j = start; j < end; j++) {
     2280          int iRow = row[j];
     2281          value += pi[iRow] * elementByColumn[j];
     2282        }
     2283      }
     2284      array[jColumn++] = value;
     2285      value = 0.0;
     2286      for (j = startNext; j < endNext; j++) {
     2287        int iRow = row[j];
     2288        value += pi[iRow] * elementByColumn[j];
     2289      }
     2290      array[jColumn] = value;
     2291    } else {
    23142292#ifdef CLP_INVESTIGATE
    2315                if (model->clpScaledMatrix())
    2316                     printf("scaledMatrix_ at %d of ClpPackedMatrix\n", __LINE__);
    2317 #endif
    2318                // scaled
    2319                const double * columnScale = model->columnScale();
    2320                int iColumn = which[0];
    2321                double value = 0.0;
    2322                double scale = columnScale[iColumn];
    2323                CoinBigIndex j;
    2324                for (j = columnStart[iColumn];
    2325                          j < columnStart[iColumn+1]; j++) {
    2326                     int iRow = row[j];
    2327                     value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
    2328                }
    2329                for (jColumn = 0; jColumn < numberToDo - 1; jColumn++) {
    2330                     int iColumn = which[jColumn+1];
    2331                     value *= scale;
    2332                     scale = columnScale[iColumn];
    2333                     CoinBigIndex start = columnStart[iColumn];
    2334                     CoinBigIndex end = columnStart[iColumn+1];
    2335                     array[jColumn] = value;
    2336                     value = 0.0;
    2337                     for (j = start; j < end; j++) {
    2338                          int iRow = row[j];
    2339                          value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
    2340                     }
    2341                }
    2342                value *= scale;
    2343                array[jColumn] = value;
    2344           }
    2345      } else if (numberToDo) {
    2346           // gaps
    2347           if (!rowScale) {
    2348                for (jColumn = 0; jColumn < numberToDo; jColumn++) {
    2349                     int iColumn = which[jColumn];
    2350                     double value = 0.0;
    2351                     CoinBigIndex j;
    2352                     for (j = columnStart[iColumn];
    2353                               j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    2354                          int iRow = row[j];
    2355                          value += pi[iRow] * elementByColumn[j];
    2356                     }
    2357                     array[jColumn] = value;
    2358                }
    2359           } else {
     2293      if (model->clpScaledMatrix())
     2294        printf("scaledMatrix_ at %d of ClpPackedMatrix\n", __LINE__);
     2295#endif
     2296      // scaled
     2297      const double *columnScale = model->columnScale();
     2298      int iColumn = which[0];
     2299      double value = 0.0;
     2300      double scale = columnScale[iColumn];
     2301      CoinBigIndex j;
     2302      for (j = columnStart[iColumn];
     2303           j < columnStart[iColumn + 1]; j++) {
     2304        int iRow = row[j];
     2305        value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
     2306      }
     2307      for (jColumn = 0; jColumn < numberToDo - 1; jColumn++) {
     2308        int iColumn = which[jColumn + 1];
     2309        value *= scale;
     2310        scale = columnScale[iColumn];
     2311        CoinBigIndex start = columnStart[iColumn];
     2312        CoinBigIndex end = columnStart[iColumn + 1];
     2313        array[jColumn] = value;
     2314        value = 0.0;
     2315        for (j = start; j < end; j++) {
     2316          int iRow = row[j];
     2317          value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
     2318        }
     2319      }
     2320      value *= scale;
     2321      array[jColumn] = value;
     2322    }
     2323  } else if (numberToDo) {
     2324    // gaps
     2325    if (!rowScale) {
     2326      for (jColumn = 0; jColumn < numberToDo; jColumn++) {
     2327        int iColumn = which[jColumn];
     2328        double value = 0.0;
     2329        CoinBigIndex j;
     2330        for (j = columnStart[iColumn];
     2331             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     2332          int iRow = row[j];
     2333          value += pi[iRow] * elementByColumn[j];
     2334        }
     2335        array[jColumn] = value;
     2336      }
     2337    } else {
    23602338#ifdef CLP_INVESTIGATE
    2361                if (model->clpScaledMatrix())
    2362                     printf("scaledMatrix_ at %d of ClpPackedMatrix - flags %d (%d) n %d\n",
    2363                            __LINE__, flags_, model->clpScaledMatrix()->flags(), numberToDo);
    2364 #endif
    2365                // scaled
    2366                const double * columnScale = model->columnScale();
    2367                for (jColumn = 0; jColumn < numberToDo; jColumn++) {
    2368                     int iColumn = which[jColumn];
    2369                     double value = 0.0;
    2370                     CoinBigIndex j;
    2371                     for (j = columnStart[iColumn];
    2372                               j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    2373                          int iRow = row[j];
    2374                          value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
    2375                     }
    2376                     value *= columnScale[iColumn];
    2377                     array[jColumn] = value;
    2378                }
    2379           }
    2380      }
     2339      if (model->clpScaledMatrix())
     2340        printf("scaledMatrix_ at %d of ClpPackedMatrix - flags %d (%d) n %d\n",
     2341          __LINE__, flags_, model->clpScaledMatrix()->flags(), numberToDo);
     2342#endif
     2343      // scaled
     2344      const double *columnScale = model->columnScale();
     2345      for (jColumn = 0; jColumn < numberToDo; jColumn++) {
     2346        int iColumn = which[jColumn];
     2347        double value = 0.0;
     2348        CoinBigIndex j;
     2349        for (j = columnStart[iColumn];
     2350             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     2351          int iRow = row[j];
     2352          value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
     2353        }
     2354        value *= columnScale[iColumn];
     2355        array[jColumn] = value;
     2356      }
     2357    }
     2358  }
    23812359}
    23822360/* Returns true if can combine transposeTimes and subsetTransposeTimes
    23832361   and if it would be faster */
    2384 bool
    2385 ClpPackedMatrix::canCombine(const ClpSimplex * model,
    2386                             const CoinIndexedVector * pi) const
    2387 {
    2388      int numberInRowArray = pi->getNumElements();
    2389      int numberRows = model->numberRows();
    2390      bool packed = pi->packedMode();
    2391      // factor should be smaller if doing both with two pi vectors
    2392      double factor = 0.30;
    2393      // We may not want to do by row if there may be cache problems
    2394      // It would be nice to find L2 cache size - for moment 512K
    2395      // Be slightly optimistic
    2396      if (numberActiveColumns_ * sizeof(double) > 1000000) {
    2397           if (numberRows * 10 < numberActiveColumns_)
    2398                factor *= 0.333333333;
    2399           else if (numberRows * 4 < numberActiveColumns_)
    2400                factor *= 0.5;
    2401           else if (numberRows * 2 < numberActiveColumns_)
    2402                factor *= 0.66666666667;
    2403           //if (model->numberIterations()%50==0)
    2404           //printf("%d nonzero\n",numberInRowArray);
    2405      }
    2406      // if not packed then bias a bit more towards by column
    2407      if (!packed)
    2408           factor *= 0.9;
    2409      // bias if columnCopy
    2410      if (columnCopy_)
    2411        factor *= 0.5;
    2412      return ((numberInRowArray > factor * numberRows || !model->rowCopy()) && !(flags_ & 2));
     2362bool ClpPackedMatrix::canCombine(const ClpSimplex *model,
     2363  const CoinIndexedVector *pi) const
     2364{
     2365  int numberInRowArray = pi->getNumElements();
     2366  int numberRows = model->numberRows();
     2367  bool packed = pi->packedMode();
     2368  // factor should be smaller if doing both with two pi vectors
     2369  double factor = 0.30;
     2370  // We may not want to do by row if there may be cache problems
     2371  // It would be nice to find L2 cache size - for moment 512K
     2372  // Be slightly optimistic
     2373  if (numberActiveColumns_ * sizeof(double) > 1000000) {
     2374    if (numberRows * 10 < numberActiveColumns_)
     2375      factor *= 0.333333333;
     2376    else if (numberRows * 4 < numberActiveColumns_)
     2377      factor *= 0.5;
     2378    else if (numberRows * 2 < numberActiveColumns_)
     2379      factor *= 0.66666666667;
     2380    //if (model->numberIterations()%50==0)
     2381    //printf("%d nonzero\n",numberInRowArray);
     2382  }
     2383  // if not packed then bias a bit more towards by column
     2384  if (!packed)
     2385    factor *= 0.9;
     2386  // bias if columnCopy
     2387  if (columnCopy_)
     2388    factor *= 0.5;
     2389  return ((numberInRowArray > factor * numberRows || !model->rowCopy()) && !(flags_ & 2));
    24132390}
    24142391#ifndef CLP_ALL_ONE_FILE
    24152392// These have to match ClpPrimalColumnSteepest version
    2416 #define reference(i)  (((reference[i>>5]>>(i&31))&1)!=0)
     2393#define reference(i) (((reference[i >> 5] >> (i & 31)) & 1) != 0)
    24172394#endif
    24182395#if ABOCA_LITE
    24192396static void
    2420 transposeTimes2UnscaledBit(clpTempInfo & info)
    2421 {
    2422   const CoinBigIndex * COIN_RESTRICT columnStart=info.start;
    2423   const int * COIN_RESTRICT row = info.row;
    2424   const double * COIN_RESTRICT elementByColumn = info.element;
    2425   double * COIN_RESTRICT array = info.infeas;
    2426   double * COIN_RESTRICT weights = const_cast<double *>(info.lower);
    2427   double * COIN_RESTRICT reducedCost = info.reducedCost;
    2428   double * COIN_RESTRICT infeas = const_cast<double *>(info.upper);
    2429   int * COIN_RESTRICT index = info.which;
    2430   unsigned int * COIN_RESTRICT reference = reinterpret_cast<unsigned int *>(info.index);
    2431   const unsigned char * COIN_RESTRICT status = info.status;
    2432   double zeroTolerance=info.tolerance;
    2433   double dualTolerance=info.dualTolerance;
    2434   const double * COIN_RESTRICT pi = info.work;
    2435   const double * COIN_RESTRICT piWeight = info.cost;
     2397transposeTimes2UnscaledBit(clpTempInfo &info)
     2398{
     2399  const CoinBigIndex *COIN_RESTRICT columnStart = info.start;
     2400  const int *COIN_RESTRICT row = info.row;
     2401  const double *COIN_RESTRICT elementByColumn = info.element;
     2402  double *COIN_RESTRICT array = info.infeas;
     2403  double *COIN_RESTRICT weights = const_cast< double * >(info.lower);
     2404  double *COIN_RESTRICT reducedCost = info.reducedCost;
     2405  double *COIN_RESTRICT infeas = const_cast< double * >(info.upper);
     2406  int *COIN_RESTRICT index = info.which;
     2407  unsigned int *COIN_RESTRICT reference = reinterpret_cast< unsigned int * >(info.index);
     2408  const unsigned char *COIN_RESTRICT status = info.status;
     2409  double zeroTolerance = info.tolerance;
     2410  double dualTolerance = info.dualTolerance;
     2411  const double *COIN_RESTRICT pi = info.work;
     2412  const double *COIN_RESTRICT piWeight = info.cost;
    24362413  double scaleFactor = info.primalRatio;
    24372414  double devex = info.changeObj;
    24382415  double referenceIn = info.theta;
    24392416  int first = info.startColumn;
    2440   int last = first+info.numberToDo;
    2441   int killDjs=info.numberInfeasibilities;
    2442   int numberNonZero=0;
    2443   double mmult[2]={1.0,-1.0};
     2417  int last = first + info.numberToDo;
     2418  int killDjs = info.numberInfeasibilities;
     2419  int numberNonZero = 0;
     2420  double mmult[2] = { 1.0, -1.0 };
    24442421  for (int iColumn = first; iColumn < last; iColumn++) {
    2445     int iStat = status[iColumn]&7;
    2446     if (iStat!=1&&iStat!=5) {
     2422    int iStat = status[iColumn] & 7;
     2423    if (iStat != 1 && iStat != 5) {
    24472424      double value = 0.0;
    24482425      CoinBigIndex start = columnStart[iColumn];
    2449       CoinBigIndex end = columnStart[iColumn+1];
    2450       int n = static_cast<int>(end - start);
    2451       const int * COIN_RESTRICT rowThis = row + start;
    2452       const double * COIN_RESTRICT elementThis = elementByColumn + start;
    2453       for (int i=0; i < n; i ++) {
    2454         int iRow = rowThis[i];
    2455         value -= pi[iRow] * elementThis[i];
     2426      CoinBigIndex end = columnStart[iColumn + 1];
     2427      int n = static_cast< int >(end - start);
     2428      const int *COIN_RESTRICT rowThis = row + start;
     2429      const double *COIN_RESTRICT elementThis = elementByColumn + start;
     2430      for (int i = 0; i < n; i++) {
     2431        int iRow = rowThis[i];
     2432        value -= pi[iRow] * elementThis[i];
    24562433      }
    24572434      if (fabs(value) > zeroTolerance) {
    2458         // and do other array
    2459         double modification = 0.0;
    2460         for (int i=0; i < n; i ++) {
    2461           int iRow = rowThis[i];
    2462           modification += piWeight[iRow] * elementThis[i];
    2463         }
    2464         double thisWeight = weights[iColumn];
    2465         double oldWeight=thisWeight;
    2466         double pivot = value * scaleFactor;
    2467         double pivotSquared = pivot * pivot;
    2468         thisWeight += pivotSquared * devex + pivot * modification;
    2469         debug3(iColumn,thisWeight,pivotSquared,devex,pivot,modification,oldWeight);
    2470         if (thisWeight < DEVEX_TRY_NORM) {
    2471           if (referenceIn < 0.0) {
    2472             // steepest
    2473             thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
    2474           } else {
    2475             // exact
    2476             thisWeight = referenceIn * pivotSquared;
    2477             if (reference(iColumn))
    2478               thisWeight += 1.0;
    2479             thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
    2480           }
    2481         }
    2482         weights[iColumn] = thisWeight;
    2483         if (!killDjs) {
    2484           value = reducedCost[iColumn]-value;
    2485           reducedCost[iColumn] = value;
    2486           if (iStat>1&&iStat<4) {
    2487             value *= mmult[iStat-2];
    2488             if (value<=dualTolerance)
    2489               value=0.0;
    2490           } else {
    2491             // free or super basic
    2492             if (fabs(value) > FREE_ACCEPT * dualTolerance) {
    2493               // we are going to bias towards free (but only if reasonable)
    2494               value *= FREE_BIAS;
    2495             } else {
    2496               value=0.0;
    2497             }
    2498           }
    2499           if (value) {
    2500             value *= value;
    2501             // store square in list
    2502             if (infeas[iColumn]) {
    2503               infeas[iColumn] = value; // already there
    2504             } else {
    2505               array[numberNonZero]=value;
    2506               index[numberNonZero++]=iColumn;
    2507             }
    2508           } else {
    2509             array[numberNonZero]=0.0;
    2510             index[numberNonZero++]=iColumn;
    2511           }
    2512         }
    2513       }
    2514     }
    2515   }
    2516   info.numberAdded=numberNonZero;
     2435        // and do other array
     2436        double modification = 0.0;
     2437        for (int i = 0; i < n; i++) {
     2438          int iRow = rowThis[i];
     2439          modification += piWeight[iRow] * elementThis[i];
     2440        }
     2441        double thisWeight = weights[iColumn];
     2442        double oldWeight = thisWeight;
     2443        double pivot = value * scaleFactor;
     2444        double pivotSquared = pivot * pivot;
     2445        thisWeight += pivotSquared * devex + pivot * modification;
     2446        debug3(iColumn, thisWeight, pivotSquared, devex, pivot, modification, oldWeight);
     2447        if (thisWeight < DEVEX_TRY_NORM) {
     2448          if (referenceIn < 0.0) {
     2449            // steepest
     2450            thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
     2451          } else {
     2452            // exact
     2453            thisWeight = referenceIn * pivotSquared;
     2454            if (reference(iColumn))
     2455              thisWeight += 1.0;
     2456            thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
     2457          }
     2458        }
     2459        weights[iColumn] = thisWeight;
     2460        if (!killDjs) {
     2461          value = reducedCost[iColumn] - value;
     2462          reducedCost[iColumn] = value;
     2463          if (iStat > 1 && iStat < 4) {
     2464            value *= mmult[iStat - 2];
     2465            if (value <= dualTolerance)
     2466              value = 0.0;
     2467          } else {
     2468            // free or super basic
     2469            if (fabs(value) > FREE_ACCEPT * dualTolerance) {
     2470              // we are going to bias towards free (but only if reasonable)
     2471              value *= FREE_BIAS;
     2472            } else {
     2473              value = 0.0;
     2474            }
     2475          }
     2476          if (value) {
     2477            value *= value;
     2478            // store square in list
     2479            if (infeas[iColumn]) {
     2480              infeas[iColumn] = value; // already there
     2481            } else {
     2482              array[numberNonZero] = value;
     2483              index[numberNonZero++] = iColumn;
     2484            }
     2485          } else {
     2486            array[numberNonZero] = 0.0;
     2487            index[numberNonZero++] = iColumn;
     2488          }
     2489        }
     2490      }
     2491    }
     2492  }
     2493  info.numberAdded = numberNonZero;
    25172494}
    25182495static void
    2519 transposeTimes2ScaledBit(clpTempInfo & info)
    2520 {
    2521   const CoinBigIndex * COIN_RESTRICT columnStart=info.start;
    2522   const int * COIN_RESTRICT row = info.row;
    2523   const double * COIN_RESTRICT elementByColumn = info.element;
    2524   double * COIN_RESTRICT array = info.infeas;
    2525   double * COIN_RESTRICT weights = const_cast<double *>(info.lower);
    2526   int * COIN_RESTRICT index = info.which;
    2527   unsigned int * COIN_RESTRICT reference = reinterpret_cast<unsigned int *>(info.index);
    2528   const unsigned char * COIN_RESTRICT status = info.status;
    2529   double zeroTolerance=info.tolerance;
    2530   double dualTolerance=info.dualTolerance;
    2531   double * COIN_RESTRICT reducedCost = info.reducedCost;
    2532   double * COIN_RESTRICT infeas = const_cast<double *>(info.upper);
    2533   const double * COIN_RESTRICT pi = info.work;
    2534   const double * COIN_RESTRICT piWeight = info.cost;
    2535   const double * COIN_RESTRICT columnScale = info.solution;
     2496transposeTimes2ScaledBit(clpTempInfo &info)
     2497{
     2498  const CoinBigIndex *COIN_RESTRICT columnStart = info.start;
     2499  const int *COIN_RESTRICT row = info.row;
     2500  const double *COIN_RESTRICT elementByColumn = info.element;
     2501  double *COIN_RESTRICT array = info.infeas;
     2502  double *COIN_RESTRICT weights = const_cast< double * >(info.lower);
     2503  int *COIN_RESTRICT index = info.which;
     2504  unsigned int *COIN_RESTRICT reference = reinterpret_cast< unsigned int * >(info.index);
     2505  const unsigned char *COIN_RESTRICT status = info.status;
     2506  double zeroTolerance = info.tolerance;
     2507  double dualTolerance = info.dualTolerance;
     2508  double *COIN_RESTRICT reducedCost = info.reducedCost;
     2509  double *COIN_RESTRICT infeas = const_cast< double * >(info.upper);
     2510  const double *COIN_RESTRICT pi = info.work;
     2511  const double *COIN_RESTRICT piWeight = info.cost;
     2512  const double *COIN_RESTRICT columnScale = info.solution;
    25362513  double scaleFactor = info.primalRatio;
    25372514  double devex = info.changeObj;
    25382515  double referenceIn = info.theta;
    25392516  int first = info.startColumn;
    2540   int last = first+info.numberToDo;
    2541   int killDjs=info.numberInfeasibilities;
    2542   int numberNonZero=0;
    2543   double mmult[2]={1.0,-1.0};
     2517  int last = first + info.numberToDo;
     2518  int killDjs = info.numberInfeasibilities;
     2519  int numberNonZero = 0;
     2520  double mmult[2] = { 1.0, -1.0 };
    25442521  for (int iColumn = first; iColumn < last; iColumn++) {
    2545     int iStat = status[iColumn]&7;
    2546     if (iStat!=1&&iStat!=5) {
     2522    int iStat = status[iColumn] & 7;
     2523    if (iStat != 1 && iStat != 5) {
    25472524      double value = 0.0;
    25482525      double scale = columnScale[iColumn];
    25492526      CoinBigIndex start = columnStart[iColumn];
    2550       CoinBigIndex end = columnStart[iColumn+1];
    2551       int n = static_cast<int>(end - start);
     2527      CoinBigIndex end = columnStart[iColumn + 1];
     2528      int n = static_cast< int >(end - start);
    25522529      bool odd = (n & 1) != 0;
    25532530      n &= ~1;
    2554       const int * COIN_RESTRICT rowThis = row + start;
    2555       const double * COIN_RESTRICT elementThis = elementByColumn + start;
    2556       for (int i=0; i < n; i += 2) {
    2557         int iRow0 = rowThis[i];
    2558         int iRow1 = rowThis[i+1];
    2559         value -= pi[iRow0] * elementThis[i];
    2560         value -= pi[iRow1] * elementThis[i+1];
     2531      const int *COIN_RESTRICT rowThis = row + start;
     2532      const double *COIN_RESTRICT elementThis = elementByColumn + start;
     2533      for (int i = 0; i < n; i += 2) {
     2534        int iRow0 = rowThis[i];
     2535        int iRow1 = rowThis[i + 1];
     2536        value -= pi[iRow0] * elementThis[i];
     2537        value -= pi[iRow1] * elementThis[i + 1];
    25612538      }
    25622539      if (odd) {
    2563         int iRow = rowThis[n];
    2564         value -= pi[iRow] * elementThis[n];
     2540        int iRow = rowThis[n];
     2541        value -= pi[iRow] * elementThis[n];
    25652542      }
    25662543      value *= scale;
    25672544      if (fabs(value) > zeroTolerance) {
    2568         // and do other array
    2569         double modification = 0.0;
    2570         for (int i=0; i < n; i += 2) {
    2571           int iRow0 = rowThis[i];
    2572           int iRow1 = rowThis[i+1];
    2573           modification += piWeight[iRow0] * elementThis[i];
    2574           modification += piWeight[iRow1] * elementThis[i+1];
    2575         }
    2576         if (odd) {
    2577           int iRow = rowThis[n];
    2578           modification += piWeight[iRow] * elementThis[n];
    2579         }
    2580         modification *= scale;
    2581         double thisWeight = weights[iColumn];
    2582         double pivot = value * scaleFactor;
    2583         double pivotSquared = pivot * pivot;
    2584         thisWeight += pivotSquared * devex + pivot * modification;
    2585         if (thisWeight < DEVEX_TRY_NORM) {
    2586           if (referenceIn < 0.0) {
    2587             // steepest
    2588             thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
    2589           } else {
    2590             // exact
    2591             thisWeight = referenceIn * pivotSquared;
    2592             if (reference(iColumn))
    2593               thisWeight += 1.0;
    2594             thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
    2595           }
    2596         }
    2597         weights[iColumn] = thisWeight;
    2598         if (!killDjs) {
    2599           value = reducedCost[iColumn]-value;
    2600           reducedCost[iColumn] = value;
    2601           if (iStat>1&&iStat<4) {
    2602             value *= mmult[iStat-2];
    2603             if (value<=dualTolerance)
    2604               value=0.0;
    2605           } else {
    2606             // free or super basic
    2607             if (fabs(value) > FREE_ACCEPT * dualTolerance) {
    2608               // we are going to bias towards free (but only if reasonable)
    2609               value *= FREE_BIAS;
    2610             } else {
    2611               value=0.0;
    2612             }
    2613           }
    2614           if (value) {
    2615             value *= value;
    2616             // store square in list
    2617             if (infeas[iColumn]) {
    2618               infeas[iColumn] = value; // already there
    2619             } else {
    2620               array[numberNonZero]=value;
    2621               index[numberNonZero++]=iColumn;
    2622             }
    2623           } else {
    2624             array[numberNonZero]=0.0;
    2625             index[numberNonZero++]=iColumn;
    2626           }
    2627         }
    2628       }
    2629     }
    2630   }
    2631   info.numberAdded=numberNonZero;
     2545        // and do other array
     2546        double modification = 0.0;
     2547        for (int i = 0; i < n; i += 2) {
     2548          int iRow0 = rowThis[i];
     2549          int iRow1 = rowThis[i + 1];
     2550          modification += piWeight[iRow0] * elementThis[i];
     2551          modification += piWeight[iRow1] * elementThis[i + 1];
     2552        }
     2553        if (odd) {
     2554          int iRow = rowThis[n];
     2555          modification += piWeight[iRow] * elementThis[n];
     2556        }
     2557        modification *= scale;
     2558        double thisWeight = weights[iColumn];
     2559        double pivot = value * scaleFactor;
     2560        double pivotSquared = pivot * pivot;
     2561        thisWeight += pivotSquared * devex + pivot * modification;
     2562        if (thisWeight < DEVEX_TRY_NORM) {
     2563          if (referenceIn < 0.0) {
     2564            // steepest
     2565            thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
     2566          } else {
     2567            // exact
     2568            thisWeight = referenceIn * pivotSquared;
     2569            if (reference(iColumn))
     2570              thisWeight += 1.0;
     2571            thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
     2572          }
     2573        }
     2574        weights[iColumn] = thisWeight;
     2575        if (!killDjs) {
     2576          value = reducedCost[iColumn] - value;
     2577          reducedCost[iColumn] = value;
     2578          if (iStat > 1 && iStat < 4) {
     2579            value *= mmult[iStat - 2];
     2580            if (value <= dualTolerance)
     2581              value = 0.0;
     2582          } else {
     2583            // free or super basic
     2584            if (fabs(value) > FREE_ACCEPT * dualTolerance) {
     2585              // we are going to bias towards free (but only if reasonable)
     2586              value *= FREE_BIAS;
     2587            } else {
     2588              value = 0.0;
     2589            }
     2590          }
     2591          if (value) {
     2592            value *= value;
     2593            // store square in list
     2594            if (infeas[iColumn]) {
     2595              infeas[iColumn] = value; // already there
     2596            } else {
     2597              array[numberNonZero] = value;
     2598              index[numberNonZero++] = iColumn;
     2599            }
     2600          } else {
     2601            array[numberNonZero] = 0.0;
     2602            index[numberNonZero++] = iColumn;
     2603          }
     2604        }
     2605      }
     2606    }
     2607  }
     2608  info.numberAdded = numberNonZero;
    26322609}
    26332610#endif
     
    26362613   new infeas in dj1 */
    26372614// Updates two arrays for steepest
    2638 int
    2639 ClpPackedMatrix::transposeTimes2(const ClpSimplex * model,
    2640                                  const CoinIndexedVector * pi1, CoinIndexedVector * dj1,
    2641                                  const CoinIndexedVector * pi2,
    2642                                  CoinIndexedVector * spare,
    2643                                  double * COIN_RESTRICT infeas,
    2644                                  double * COIN_RESTRICT reducedCost,
    2645                                  double referenceIn, double devex,
    2646                                  // Array for exact devex to say what is in reference framework
    2647                                  unsigned int * COIN_RESTRICT reference,
    2648                                  double * COIN_RESTRICT weights, double scaleFactor)
    2649 {
    2650      int returnCode=0;
    2651      // put row of tableau in dj1
    2652      double * COIN_RESTRICT pi = pi1->denseVector();
    2653      int numberNonZero = 0;
    2654      int * COIN_RESTRICT index = dj1->getIndices();
    2655      double * COIN_RESTRICT array = dj1->denseVector();
    2656      int numberInRowArray = pi1->getNumElements();
    2657      double zeroTolerance = model->zeroTolerance();
    2658      double dualTolerance = model->currentDualTolerance();
    2659      // we can't really trust infeasibilities if there is dual error
    2660      // this coding has to mimic coding in checkDualSolution
    2661      double error = CoinMin(1.0e-2, model->largestDualError());
    2662      // allow tolerance at least slightly bigger than standard
    2663      dualTolerance = dualTolerance +  error;
    2664      bool packed = pi1->packedMode();
    2665      // do by column
    2666      int iColumn;
    2667      // get matrix data pointers
    2668      const int * COIN_RESTRICT row = matrix_->getIndices();
    2669      const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
    2670      const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
    2671      const double * COIN_RESTRICT rowScale = model->rowScale();
    2672      assert (!spare->getNumElements());
    2673      assert (numberActiveColumns_ > 0);
    2674      double * COIN_RESTRICT piWeight = pi2->denseVector();
    2675      assert (!pi2->packedMode());
    2676      bool killDjs = (scaleFactor == 0.0);
    2677      if (!scaleFactor)
    2678           scaleFactor = 1.0;
    2679      if (packed) {
    2680           // need to expand pi into y
    2681           assert(spare->capacity() >= model->numberRows());
    2682           double * COIN_RESTRICT piOld = pi;
    2683           pi = spare->denseVector();
    2684           const int * COIN_RESTRICT whichRow = pi1->getIndices();
    2685           int i;
    2686           ClpPackedMatrix * scaledMatrix = model->clpScaledMatrix();
    2687           if (rowScale && scaledMatrix) {
    2688                rowScale = NULL;
    2689                // get matrix data pointers
    2690                row = scaledMatrix->getIndices();
    2691                columnStart = scaledMatrix->getVectorStarts();
    2692                elementByColumn = scaledMatrix->getElements();
    2693           }
    2694           if (!rowScale) {
    2695                // modify pi so can collapse to one loop
    2696                for (i = 0; i < numberInRowArray; i++) {
    2697                     int iRow = whichRow[i];
    2698                     pi[iRow] = piOld[i];
    2699                }
    2700                if (!columnCopy_ || killDjs) {
    2701                  if (infeas)
    2702                    returnCode=1;
     2615int ClpPackedMatrix::transposeTimes2(const ClpSimplex *model,
     2616  const CoinIndexedVector *pi1, CoinIndexedVector *dj1,
     2617  const CoinIndexedVector *pi2,
     2618  CoinIndexedVector *spare,
     2619  double *COIN_RESTRICT infeas,
     2620  double *COIN_RESTRICT reducedCost,
     2621  double referenceIn, double devex,
     2622  // Array for exact devex to say what is in reference framework
     2623  unsigned int *COIN_RESTRICT reference,
     2624  double *COIN_RESTRICT weights, double scaleFactor)
     2625{
     2626  int returnCode = 0;
     2627  // put row of tableau in dj1
     2628  double *COIN_RESTRICT pi = pi1->denseVector();
     2629  int numberNonZero = 0;
     2630  int *COIN_RESTRICT index = dj1->getIndices();
     2631  double *COIN_RESTRICT array = dj1->denseVector();
     2632  int numberInRowArray = pi1->getNumElements();
     2633  double zeroTolerance = model->zeroTolerance();
     2634  double dualTolerance = model->currentDualTolerance();
     2635  // we can't really trust infeasibilities if there is dual error
     2636  // this coding has to mimic coding in checkDualSolution
     2637  double error = CoinMin(1.0e-2, model->largestDualError());
     2638  // allow tolerance at least slightly bigger than standard
     2639  dualTolerance = dualTolerance + error;
     2640  bool packed = pi1->packedMode();
     2641  // do by column
     2642  int iColumn;
     2643  // get matrix data pointers
     2644  const int *COIN_RESTRICT row = matrix_->getIndices();
     2645  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     2646  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
     2647  const double *COIN_RESTRICT rowScale = model->rowScale();
     2648  assert(!spare->getNumElements());
     2649  assert(numberActiveColumns_ > 0);
     2650  double *COIN_RESTRICT piWeight = pi2->denseVector();
     2651  assert(!pi2->packedMode());
     2652  bool killDjs = (scaleFactor == 0.0);
     2653  if (!scaleFactor)
     2654    scaleFactor = 1.0;
     2655  if (packed) {
     2656    // need to expand pi into y
     2657    assert(spare->capacity() >= model->numberRows());
     2658    double *COIN_RESTRICT piOld = pi;
     2659    pi = spare->denseVector();
     2660    const int *COIN_RESTRICT whichRow = pi1->getIndices();
     2661    int i;
     2662    ClpPackedMatrix *scaledMatrix = model->clpScaledMatrix();
     2663    if (rowScale && scaledMatrix) {
     2664      rowScale = NULL;
     2665      // get matrix data pointers
     2666      row = scaledMatrix->getIndices();
     2667      columnStart = scaledMatrix->getVectorStarts();
     2668      elementByColumn = scaledMatrix->getElements();
     2669    }
     2670    if (!rowScale) {
     2671      // modify pi so can collapse to one loop
     2672      for (i = 0; i < numberInRowArray; i++) {
     2673        int iRow = whichRow[i];
     2674        pi[iRow] = piOld[i];
     2675      }
     2676      if (!columnCopy_ || killDjs) {
     2677        if (infeas)
     2678          returnCode = 1;
    27032679#if ABOCA_LITE
    2704                int numberThreads=abcState();
    2705                if (numberThreads) {
    2706                  clpTempInfo info[ABOCA_LITE];
    2707                  int chunk = (numberActiveColumns_+numberThreads-1)/numberThreads;
    2708                  int n=0;
    2709                  for (int i=0;i<numberThreads;i++) {
    2710                    info[i].theta = referenceIn;
    2711                    info[i].changeObj = devex;
    2712                    info[i].primalRatio = scaleFactor;
    2713                    info[i].lower = weights;
    2714                    info[i].reducedCost = reducedCost;
    2715                    info[i].upper = infeas;
    2716                    info[i].which = index+n;
    2717                    info[i].infeas = array+n;
    2718                    info[i].index = reinterpret_cast<int *>(reference);
    2719                    info[i].work=const_cast<double *>(pi);
    2720                    info[i].cost=const_cast<double *>(piWeight);
    2721                    info[i].status=const_cast<unsigned char *>(model->statusArray());
    2722                    info[i].startColumn=n;
    2723                    info[i].element=elementByColumn;
    2724                    info[i].start=columnStart;
    2725                    info[i].row=row;
    2726                    info[i].numberToDo=CoinMin(chunk,numberActiveColumns_-n);
    2727                    info[i].tolerance=zeroTolerance;
    2728                    info[i].dualTolerance=dualTolerance;
    2729                    info[i].numberInfeasibilities=killDjs ? 1 : 0;
    2730                    n += chunk;
    2731                  }
    2732                  for (int i=0;i<numberThreads;i++) {
    2733                    cilk_spawn transposeTimes2UnscaledBit(info[i]);
    2734                  }
    2735                  cilk_sync;
    2736                  for (int i=0;i<numberThreads;i++) {
    2737                    numberNonZero += info[i].numberAdded;
    2738                  }
    2739                  moveAndZero(info,2,NULL);
    2740                } else {
    2741 #endif
    2742                     CoinBigIndex j;
    2743                     CoinBigIndex end = columnStart[0];
    2744                     for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
    2745                       CoinBigIndex start = end;
    2746                       end = columnStart[iColumn+1];
    2747                       ClpSimplex::Status status = model->getStatus(iColumn);
    2748                       if (status == ClpSimplex::basic || status == ClpSimplex::isFixed) continue;
    2749                       double value = 0.0;
    2750                       for (j = start; j < end; j++) {
    2751                         int iRow = row[j];
    2752                         value -= pi[iRow] * elementByColumn[j];
    2753                       }
    2754                       if (fabs(value) > zeroTolerance) {
    2755                         // and do other array
    2756                         double modification = 0.0;
    2757                         for (j = start; j < end; j++) {
    2758                           int iRow = row[j];
    2759                           modification += piWeight[iRow] * elementByColumn[j];
    2760                         }
    2761                         double thisWeight = weights[iColumn];
    2762                         double oldWeight=thisWeight;
    2763                         double pivot = value * scaleFactor;
    2764                         double pivotSquared = pivot * pivot;
    2765                         thisWeight += pivotSquared * devex + pivot * modification;
    2766                         debug3(iColumn,thisWeight,pivotSquared,devex,pivot,modification,oldWeight);
    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 (reference(iColumn))
    2775                               thisWeight += 1.0;
    2776                             thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
    2777                           }
    2778                         }
    2779                         weights[iColumn] = thisWeight;
    2780                         if (!killDjs) {
    2781                           value = reducedCost[iColumn]-value;
    2782                           reducedCost[iColumn] = value;
    2783                           // simplify status
    2784                           ClpSimplex::Status status = model->getStatus(iColumn);
    2785                          
    2786                           switch(status) {
    2787                            
    2788                           case ClpSimplex::basic:
    2789                           case ClpSimplex::isFixed:
    2790                             break;
    2791                           case ClpSimplex::isFree:
    2792                           case ClpSimplex::superBasic:
    2793                             if (fabs(value) > FREE_ACCEPT * dualTolerance) {
    2794                               // we are going to bias towards free (but only if reasonable)
    2795                               value *= FREE_BIAS;
    2796                               value *= value;
    2797                               // store square in list
    2798                               if (infeas[iColumn]) {
    2799                                 infeas[iColumn] = value; // already there
    2800                               } else {
    2801                                 array[numberNonZero]=value;
    2802                                 index[numberNonZero++]=iColumn;
    2803                               }
    2804                             } else {
    2805                               array[numberNonZero]=0.0;
    2806                               index[numberNonZero++]=iColumn;
    2807                             }
    2808                             break;
    2809                           case ClpSimplex::atUpperBound:
    2810                             if (value > dualTolerance) {
    2811                               value *= value;
    2812                               // store square in list
    2813                               if (infeas[iColumn]) {
    2814                                 infeas[iColumn] = value; // already there
    2815                               } else {
    2816                                 array[numberNonZero]=value;
    2817                                 index[numberNonZero++]=iColumn;
    2818                               }
    2819                             } else {
    2820                               array[numberNonZero]=0.0;
    2821                               index[numberNonZero++]=iColumn;
    2822                             }
    2823                             break;
    2824                           case ClpSimplex::atLowerBound:
    2825                             if (value < -dualTolerance) {
    2826                               value *= value;
    2827                               // store square in list
    2828                               if (infeas[iColumn]) {
    2829                                 infeas[iColumn] = value; // already there
    2830                               } else {
    2831                                 array[numberNonZero]=value;
    2832                                 index[numberNonZero++]=iColumn;
    2833                               }
    2834                             } else {
    2835                               array[numberNonZero]=0.0;
    2836                               index[numberNonZero++]=iColumn;
    2837                             }
    2838                           }
    2839                         }
    2840                       }
     2680        int numberThreads = abcState();
     2681        if (numberThreads) {
     2682          clpTempInfo info[ABOCA_LITE];
     2683          int chunk = (numberActiveColumns_ + numberThreads - 1) / numberThreads;
     2684          int n = 0;
     2685          for (int i = 0; i < numberThreads; i++) {
     2686            info[i].theta = referenceIn;
     2687            info[i].changeObj = devex;
     2688            info[i].primalRatio = scaleFactor;
     2689            info[i].lower = weights;
     2690            info[i].reducedCost = reducedCost;
     2691            info[i].upper = infeas;
     2692            info[i].which = index + n;
     2693            info[i].infeas = array + n;
     2694            info[i].index = reinterpret_cast< int * >(reference);
     2695            info[i].work = const_cast< double * >(pi);
     2696            info[i].cost = const_cast< double * >(piWeight);
     2697            info[i].status = const_cast< unsigned char * >(model->statusArray());
     2698            info[i].startColumn = n;
     2699            info[i].element = elementByColumn;
     2700            info[i].start = columnStart;
     2701            info[i].row = row;
     2702            info[i].numberToDo = CoinMin(chunk, numberActiveColumns_ - n);
     2703            info[i].tolerance = zeroTolerance;
     2704            info[i].dualTolerance = dualTolerance;
     2705            info[i].numberInfeasibilities = killDjs ? 1 : 0;
     2706            n += chunk;
     2707          }
     2708          for (int i = 0; i < numberThreads; i++) {
     2709            cilk_spawn transposeTimes2UnscaledBit(info[i]);
     2710          }
     2711          cilk_sync;
     2712          for (int i = 0; i < numberThreads; i++) {
     2713            numberNonZero += info[i].numberAdded;
     2714          }
     2715          moveAndZero(info, 2, NULL);
     2716        } else {
     2717#endif
     2718          CoinBigIndex j;
     2719          CoinBigIndex end = columnStart[0];
     2720          for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
     2721            CoinBigIndex start = end;
     2722            end = columnStart[iColumn + 1];
     2723            ClpSimplex::Status status = model->getStatus(iColumn);
     2724            if (status == ClpSimplex::basic || status == ClpSimplex::isFixed)
     2725              continue;
     2726            double value = 0.0;
     2727            for (j = start; j < end; j++) {
     2728              int iRow = row[j];
     2729              value -= pi[iRow] * elementByColumn[j];
     2730            }
     2731            if (fabs(value) > zeroTolerance) {
     2732              // and do other array
     2733              double modification = 0.0;
     2734              for (j = start; j < end; j++) {
     2735                int iRow = row[j];
     2736                modification += piWeight[iRow] * elementByColumn[j];
     2737              }
     2738              double thisWeight = weights[iColumn];
     2739              double oldWeight = thisWeight;
     2740              double pivot = value * scaleFactor;
     2741              double pivotSquared = pivot * pivot;
     2742              thisWeight += pivotSquared * devex + pivot * modification;
     2743              debug3(iColumn, thisWeight, pivotSquared, devex, pivot, modification, oldWeight);
     2744              if (thisWeight < DEVEX_TRY_NORM) {
     2745                if (referenceIn < 0.0) {
     2746                  // steepest
     2747                  thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
     2748                } else {
     2749                  // exact
     2750                  thisWeight = referenceIn * pivotSquared;
     2751                  if (reference(iColumn))
     2752                    thisWeight += 1.0;
     2753                  thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
     2754                }
     2755              }
     2756              weights[iColumn] = thisWeight;
     2757              if (!killDjs) {
     2758                value = reducedCost[iColumn] - value;
     2759                reducedCost[iColumn] = value;
     2760                // simplify status
     2761                ClpSimplex::Status status = model->getStatus(iColumn);
     2762
     2763                switch (status) {
     2764
     2765                case ClpSimplex::basic:
     2766                case ClpSimplex::isFixed:
     2767                  break;
     2768                case ClpSimplex::isFree:
     2769                case ClpSimplex::superBasic:
     2770                  if (fabs(value) > FREE_ACCEPT * dualTolerance) {
     2771                    // we are going to bias towards free (but only if reasonable)
     2772                    value *= FREE_BIAS;
     2773                    value *= value;
     2774                    // store square in list
     2775                    if (infeas[iColumn]) {
     2776                      infeas[iColumn] = value; // already there
     2777                    } else {
     2778                      array[numberNonZero] = value;
     2779                      index[numberNonZero++] = iColumn;
    28412780                    }
     2781                  } else {
     2782                    array[numberNonZero] = 0.0;
     2783                    index[numberNonZero++] = iColumn;
     2784                  }
     2785                  break;
     2786                case ClpSimplex::atUpperBound:
     2787                  if (value > dualTolerance) {
     2788                    value *= value;
     2789                    // store square in list
     2790                    if (infeas[iColumn]) {
     2791                      infeas[iColumn] = value; // already there
     2792                    } else {
     2793                      array[numberNonZero] = value;
     2794                      index[numberNonZero++] = iColumn;
     2795                    }
     2796                  } else {
     2797                    array[numberNonZero] = 0.0;
     2798                    index[numberNonZero++] = iColumn;
     2799                  }
     2800                  break;
     2801                case ClpSimplex::atLowerBound:
     2802                  if (value < -dualTolerance) {
     2803                    value *= value;
     2804                    // store square in list
     2805                    if (infeas[iColumn]) {
     2806                      infeas[iColumn] = value; // already there
     2807                    } else {
     2808                      array[numberNonZero] = value;
     2809                      index[numberNonZero++] = iColumn;
     2810                    }
     2811                  } else {
     2812                    array[numberNonZero] = 0.0;
     2813                    index[numberNonZero++] = iColumn;
     2814                  }
     2815                }
     2816              }
     2817            }
     2818          }
    28422819#if ABOCA_LITE
    2843             }
    2844 #endif
    2845                } else {
    2846                     // use special column copy
    2847                     // reset back
    2848                     assert (!killDjs);
    2849                     //if (killDjs)
    2850                     //   scaleFactor = 0.0;
    2851                 if (infeas)
    2852                    returnCode=1;
    2853                     columnCopy_->transposeTimes2(model, pi, dj1, piWeight,
    2854                                                 infeas, reducedCost,
    2855                                                 referenceIn, devex,
    2856                                                  reference, weights, scaleFactor);
    2857                     numberNonZero = dj1->getNumElements();
    2858                }
    2859           } else {
    2860                // scaled
    2861                // modify pi so can collapse to one loop
    2862                for (i = 0; i < numberInRowArray; i++) {
    2863                     int iRow = whichRow[i];
    2864                     pi[iRow] = piOld[i] * rowScale[iRow];
    2865                }
    2866                // can also scale piWeight as not used again
    2867                int numberWeight = pi2->getNumElements();
    2868                const int * indexWeight = pi2->getIndices();
    2869                for (i = 0; i < numberWeight; i++) {
    2870                     int iRow = indexWeight[i];
    2871                     piWeight[iRow] *= rowScale[iRow];
    2872                }
    2873                if (!columnCopy_ || killDjs) {
    2874                 if (infeas)
    2875                    returnCode=1;
    2876                  const double * COIN_RESTRICT columnScale = model->columnScale();
     2820        }
     2821#endif
     2822      } else {
     2823        // use special column copy
     2824        // reset back
     2825        assert(!killDjs);
     2826        //if (killDjs)
     2827        //   scaleFactor = 0.0;
     2828        if (infeas)
     2829          returnCode = 1;
     2830        columnCopy_->transposeTimes2(model, pi, dj1, piWeight,
     2831          infeas, reducedCost,
     2832          referenceIn, devex,
     2833          reference, weights, scaleFactor);
     2834        numberNonZero = dj1->getNumElements();
     2835      }
     2836    } else {
     2837      // scaled
     2838      // modify pi so can collapse to one loop
     2839      for (i = 0; i < numberInRowArray; i++) {
     2840        int iRow = whichRow[i];
     2841        pi[iRow] = piOld[i] * rowScale[iRow];
     2842      }
     2843      // can also scale piWeight as not used again
     2844      int numberWeight = pi2->getNumElements();
     2845      const int *indexWeight = pi2->getIndices();
     2846      for (i = 0; i < numberWeight; i++) {
     2847        int iRow = indexWeight[i];
     2848        piWeight[iRow] *= rowScale[iRow];
     2849      }
     2850      if (!columnCopy_ || killDjs) {
     2851        if (infeas)
     2852          returnCode = 1;
     2853        const double *COIN_RESTRICT columnScale = model->columnScale();
    28772854#if ABOCA_LITE
    2878                  int numberThreads=abcState();
    2879                  if (numberThreads) {
    2880                    clpTempInfo info[ABOCA_LITE];
    2881                    int chunk = (numberActiveColumns_+numberThreads-1)/numberThreads;
    2882                    int n=0;
    2883                    for (int i=0;i<numberThreads;i++) {
    2884                      info[i].theta = referenceIn;
    2885                      info[i].changeObj = devex;
    2886                      info[i].primalRatio = scaleFactor;
    2887                      info[i].lower = weights;
    2888                      info[i].reducedCost = reducedCost;
    2889                      info[i].upper = infeas;
    2890                      info[i].solution = const_cast<double *>(columnScale);
    2891                      info[i].which = index+n;
    2892                      info[i].infeas = array+n;
    2893                      info[i].index = reinterpret_cast<int *>(reference);
    2894                      info[i].work=const_cast<double *>(pi);
    2895                      info[i].cost=const_cast<double *>(piWeight);
    2896                      info[i].status=const_cast<unsigned char *>(model->statusArray());
    2897                      info[i].startColumn=n;
    2898                      info[i].element=elementByColumn;
    2899                      info[i].start=columnStart;
    2900                      info[i].row=row;
    2901                      info[i].numberToDo=CoinMin(chunk,numberActiveColumns_-n);
    2902                      info[i].tolerance=zeroTolerance;
    2903                      info[i].dualTolerance=dualTolerance;
    2904                      info[i].numberInfeasibilities=killDjs ? 1 : 0;
    2905                      n += chunk;
    2906                    }
    2907                    for (int i=0;i<numberThreads;i++) {
    2908                      cilk_spawn transposeTimes2ScaledBit(info[i]);
    2909                    }
    2910                    cilk_sync;
    2911                    for (int i=0;i<numberThreads;i++) {
    2912                      numberNonZero += info[i].numberAdded;
    2913                    }
    2914                    moveAndZero(info,2,NULL);
    2915                    if (infeas) {
    2916                      returnCode=1;
    2917                      dj1->setNumElements(numberNonZero);
    2918                    }
    2919                  } else {
    2920 #endif
    2921                     CoinBigIndex j;
    2922                     CoinBigIndex end = columnStart[0];
    2923                     for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
    2924                          CoinBigIndex start = end;
    2925                          end = columnStart[iColumn+1];
    2926                          ClpSimplex::Status status = model->getStatus(iColumn);
    2927                          if (status == ClpSimplex::basic || status == ClpSimplex::isFixed) continue;
    2928                          double scale = columnScale[iColumn];
    2929                          double value = 0.0;
    2930                          for (j = start; j < end; j++) {
    2931                               int iRow = row[j];
    2932                               value -= pi[iRow] * elementByColumn[j];
    2933                          }
    2934                          value *= scale;
    2935                          if (fabs(value) > zeroTolerance) {
    2936                               double modification = 0.0;
    2937                               for (j = start; j < end; j++) {
    2938                                    int iRow = row[j];
    2939                                    modification += piWeight[iRow] * elementByColumn[j];
    2940                               }
    2941                               modification *= scale;
    2942                               double thisWeight = weights[iColumn];
    2943                               double pivot = value * scaleFactor;
    2944                               double pivotSquared = pivot * pivot;
    2945                               thisWeight += pivotSquared * devex + pivot * modification;
    2946                               if (thisWeight < DEVEX_TRY_NORM) {
    2947                                    if (referenceIn < 0.0) {
    2948                                         // steepest
    2949                                         thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
    2950                                    } else {
    2951                                         // exact
    2952                                         thisWeight = referenceIn * pivotSquared;
    2953                                         if (reference(iColumn))
    2954                                              thisWeight += 1.0;
    2955                                         thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
    2956                                    }
    2957                               }
    2958                               weights[iColumn] = thisWeight;
    2959                               if (!killDjs) {
    2960                                 value = reducedCost[iColumn]-value;
    2961                                 reducedCost[iColumn] = value;
    2962                                 // simplify status
    2963                                 ClpSimplex::Status status = model->getStatus(iColumn);
    2964                                
    2965                                 switch(status) {
    2966                                  
    2967                                 case ClpSimplex::basic:
    2968                                 case ClpSimplex::isFixed:
    2969                                   break;
    2970                                 case ClpSimplex::isFree:
    2971                                 case ClpSimplex::superBasic:
    2972                                   if (fabs(value) > FREE_ACCEPT * dualTolerance) {
    2973                                     // we are going to bias towards free (but only if reasonable)
    2974                                     value *= FREE_BIAS;
    2975                                     value *= value;
    2976                                     // store square in list
    2977                                     if (infeas[iColumn]) {
    2978                                       infeas[iColumn] = value; // already there
    2979                                     } else {
    2980                                       array[numberNonZero]=value;
    2981                                       index[numberNonZero++]=iColumn;
    2982                                     }
    2983                                   } else {
    2984                                     array[numberNonZero]=0.0;
    2985                                     index[numberNonZero++]=iColumn;
    2986                                   }
    2987                                   break;
    2988                                 case ClpSimplex::atUpperBound:
    2989                                   if (value > dualTolerance) {
    2990                                     value *= value;
    2991                                     // store square in list
    2992                                     if (infeas[iColumn]) {
    2993                                       infeas[iColumn] = value; // already there
    2994                                     } else {
    2995                                       array[numberNonZero]=value;
    2996                                       index[numberNonZero++]=iColumn;
    2997                                     }
    2998                                   } else {
    2999                                     array[numberNonZero]=0.0;
    3000                                     index[numberNonZero++]=iColumn;
    3001                                   }
    3002                                   break;
    3003                                 case ClpSimplex::atLowerBound:
    3004                                   if (value < -dualTolerance) {
    3005                                     value *= value;
    3006                                     // store square in list
    3007                                     if (infeas[iColumn]) {
    3008                                       infeas[iColumn] = value; // already there
    3009                                     } else {
    3010                                       array[numberNonZero]=value;
    3011                                       index[numberNonZero++]=iColumn;
    3012                                     }
    3013                                   } else {
    3014                                     array[numberNonZero]=0.0;
    3015                                     index[numberNonZero++]=iColumn;
    3016                                   }
    3017                                 }
    3018                               }
    3019                          }
     2855        int numberThreads = abcState();
     2856        if (numberThreads) {
     2857          clpTempInfo info[ABOCA_LITE];
     2858          int chunk = (numberActiveColumns_ + numberThreads - 1) / numberThreads;
     2859          int n = 0;
     2860          for (int i = 0; i < numberThreads; i++) {
     2861            info[i].theta = referenceIn;
     2862            info[i].changeObj = devex;
     2863            info[i].primalRatio = scaleFactor;
     2864            info[i].lower = weights;
     2865            info[i].reducedCost = reducedCost;
     2866            info[i].upper = infeas;
     2867            info[i].solution = const_cast< double * >(columnScale);
     2868            info[i].which = index + n;
     2869            info[i].infeas = array + n;
     2870            info[i].index = reinterpret_cast< int * >(reference);
     2871            info[i].work = const_cast< double * >(pi);
     2872            info[i].cost = const_cast< double * >(piWeight);
     2873            info[i].status = const_cast< unsigned char * >(model->statusArray());
     2874            info[i].startColumn = n;
     2875            info[i].element = elementByColumn;
     2876            info[i].start = columnStart;
     2877            info[i].row = row;
     2878            info[i].numberToDo = CoinMin(chunk, numberActiveColumns_ - n);
     2879            info[i].tolerance = zeroTolerance;
     2880            info[i].dualTolerance = dualTolerance;
     2881            info[i].numberInfeasibilities = killDjs ? 1 : 0;
     2882            n += chunk;
     2883          }
     2884          for (int i = 0; i < numberThreads; i++) {
     2885            cilk_spawn transposeTimes2ScaledBit(info[i]);
     2886          }
     2887          cilk_sync;
     2888          for (int i = 0; i < numberThreads; i++) {
     2889            numberNonZero += info[i].numberAdded;
     2890          }
     2891          moveAndZero(info, 2, NULL);
     2892          if (infeas) {
     2893            returnCode = 1;
     2894            dj1->setNumElements(numberNonZero);
     2895          }
     2896        } else {
     2897#endif
     2898          CoinBigIndex j;
     2899          CoinBigIndex end = columnStart[0];
     2900          for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
     2901            CoinBigIndex start = end;
     2902            end = columnStart[iColumn + 1];
     2903            ClpSimplex::Status status = model->getStatus(iColumn);
     2904            if (status == ClpSimplex::basic || status == ClpSimplex::isFixed)
     2905              continue;
     2906            double scale = columnScale[iColumn];
     2907            double value = 0.0;
     2908            for (j = start; j < end; j++) {
     2909              int iRow = row[j];
     2910              value -= pi[iRow] * elementByColumn[j];
     2911            }
     2912            value *= scale;
     2913            if (fabs(value) > zeroTolerance) {
     2914              double modification = 0.0;
     2915              for (j = start; j < end; j++) {
     2916                int iRow = row[j];
     2917                modification += piWeight[iRow] * elementByColumn[j];
     2918              }
     2919              modification *= scale;
     2920              double thisWeight = weights[iColumn];
     2921              double pivot = value * scaleFactor;
     2922              double pivotSquared = pivot * pivot;
     2923              thisWeight += pivotSquared * devex + pivot * modification;
     2924              if (thisWeight < DEVEX_TRY_NORM) {
     2925                if (referenceIn < 0.0) {
     2926                  // steepest
     2927                  thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
     2928                } else {
     2929                  // exact
     2930                  thisWeight = referenceIn * pivotSquared;
     2931                  if (reference(iColumn))
     2932                    thisWeight += 1.0;
     2933                  thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
     2934                }
     2935              }
     2936              weights[iColumn] = thisWeight;
     2937              if (!killDjs) {
     2938                value = reducedCost[iColumn] - value;
     2939                reducedCost[iColumn] = value;