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

formatting

File:
1 edited

Legend:

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

    r2271 r2385  
    33// Corporation and others.  All Rights Reserved.
    44// This code is licensed under the terms of the Eclipse Public License (EPL).
    5 
    65
    76#include <cstdio>
     
    2827// Default Constructor
    2928//-------------------------------------------------------------------
    30 ClpNetworkMatrix::ClpNetworkMatrix ()
    31      : ClpMatrixBase()
    32 {
    33      setType(11);
    34      matrix_ = NULL;
    35      lengths_ = NULL;
    36      indices_ = NULL;
    37      numberRows_ = 0;
    38      numberColumns_ = 0;
    39      trueNetwork_ = false;
     29ClpNetworkMatrix::ClpNetworkMatrix()
     30  : ClpMatrixBase()
     31{
     32  setType(11);
     33  matrix_ = NULL;
     34  lengths_ = NULL;
     35  indices_ = NULL;
     36  numberRows_ = 0;
     37  numberColumns_ = 0;
     38  trueNetwork_ = false;
    4039}
    4140
    4241/* Constructor from two arrays */
    43 ClpNetworkMatrix::ClpNetworkMatrix(int numberColumns, const int * head,
    44                                    const int * tail)
    45      : ClpMatrixBase()
    46 {
    47      setType(11);
    48      matrix_ = NULL;
    49      lengths_ = NULL;
    50      indices_ = new int[2*numberColumns];;
    51      numberRows_ = -1;
    52      numberColumns_ = numberColumns;
    53      trueNetwork_ = true;
    54      int iColumn;
    55      CoinBigIndex j = 0;
    56      for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
    57           int iRow = head[iColumn];
    58           numberRows_ = CoinMax(numberRows_, iRow);
    59           indices_[j] = iRow;
    60           iRow = tail[iColumn];
    61           numberRows_ = CoinMax(numberRows_, iRow);
    62           indices_[j+1] = iRow;
    63      }
    64      numberRows_++;
     42ClpNetworkMatrix::ClpNetworkMatrix(int numberColumns, const int *head,
     43  const int *tail)
     44  : ClpMatrixBase()
     45{
     46  setType(11);
     47  matrix_ = NULL;
     48  lengths_ = NULL;
     49  indices_ = new int[2 * numberColumns];
     50  ;
     51  numberRows_ = -1;
     52  numberColumns_ = numberColumns;
     53  trueNetwork_ = true;
     54  int iColumn;
     55  CoinBigIndex j = 0;
     56  for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
     57    int iRow = head[iColumn];
     58    numberRows_ = CoinMax(numberRows_, iRow);
     59    indices_[j] = iRow;
     60    iRow = tail[iColumn];
     61    numberRows_ = CoinMax(numberRows_, iRow);
     62    indices_[j + 1] = iRow;
     63  }
     64  numberRows_++;
    6565}
    6666//-------------------------------------------------------------------
    6767// Copy constructor
    6868//-------------------------------------------------------------------
    69 ClpNetworkMatrix::ClpNetworkMatrix (const ClpNetworkMatrix & rhs)
    70      : ClpMatrixBase(rhs)
    71 {
    72      matrix_ = NULL;
    73      lengths_ = NULL;
    74      indices_ = NULL;
    75      numberRows_ = rhs.numberRows_;
    76      numberColumns_ = rhs.numberColumns_;
    77      trueNetwork_ = rhs.trueNetwork_;
    78      if (numberColumns_) {
    79           indices_ = new int [ 2*numberColumns_];
    80           CoinMemcpyN(rhs.indices_, 2 * numberColumns_, indices_);
    81      }
    82      int numberRows = getNumRows();
    83      if (rhs.rhsOffset_ && numberRows) {
    84           rhsOffset_ = ClpCopyOfArray(rhs.rhsOffset_, numberRows);
    85      } else {
    86           rhsOffset_ = NULL;
    87      }
    88 }
    89 
    90 ClpNetworkMatrix::ClpNetworkMatrix (const CoinPackedMatrix & rhs)
    91      : ClpMatrixBase()
    92 {
    93      setType(11);
    94      matrix_ = NULL;
    95      lengths_ = NULL;
    96      indices_ = NULL;
    97      int iColumn;
    98      assert (rhs.isColOrdered());
    99      // get matrix data pointers
    100      const int * row = rhs.getIndices();
    101      const CoinBigIndex * columnStart = rhs.getVectorStarts();
    102      const int * columnLength = rhs.getVectorLengths();
    103      const double * elementByColumn = rhs.getElements();
    104      numberColumns_ = rhs.getNumCols();
    105      int goodNetwork = 1;
    106      numberRows_ = -1;
    107      indices_ = new int[2*numberColumns_];
    108      CoinBigIndex j = 0;
    109      for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
    110           CoinBigIndex k = columnStart[iColumn];
    111           int iRow;
    112           switch (columnLength[iColumn]) {
    113           case 0:
    114                goodNetwork = -1; // not classic network
    115                indices_[j] = -1;
    116                indices_[j+1] = -1;
    117                break;
    118 
    119           case 1:
    120                goodNetwork = -1; // not classic network
    121                if (fabs(elementByColumn[k] - 1.0) < 1.0e-10) {
    122                     indices_[j] = -1;
    123                     iRow = row[k];
    124                     numberRows_ = CoinMax(numberRows_, iRow);
    125                     indices_[j+1] = iRow;
    126                } else if (fabs(elementByColumn[k] + 1.0) < 1.0e-10) {
    127                     indices_[j+1] = -1;
    128                     iRow = row[k];
    129                     numberRows_ = CoinMax(numberRows_, iRow);
    130                     indices_[j] = iRow;
    131                } else {
    132                     goodNetwork = 0; // not a network
    133                }
    134                break;
    135 
    136           case 2:
    137                if (fabs(elementByColumn[k] - 1.0) < 1.0e-10) {
    138                     if (fabs(elementByColumn[k+1] + 1.0) < 1.0e-10) {
    139                          iRow = row[k];
    140                          numberRows_ = CoinMax(numberRows_, iRow);
    141                          indices_[j+1] = iRow;
    142                          iRow = row[k+1];
    143                          numberRows_ = CoinMax(numberRows_, iRow);
    144                          indices_[j] = iRow;
    145                     } else {
    146                          goodNetwork = 0; // not a network
    147                     }
    148                } else if (fabs(elementByColumn[k] + 1.0) < 1.0e-10) {
    149                     if (fabs(elementByColumn[k+1] - 1.0) < 1.0e-10) {
    150                          iRow = row[k];
    151                          numberRows_ = CoinMax(numberRows_, iRow);
    152                          indices_[j] = iRow;
    153                          iRow = row[k+1];
    154                          numberRows_ = CoinMax(numberRows_, iRow);
    155                          indices_[j+1] = iRow;
    156                     } else {
    157                          goodNetwork = 0; // not a network
    158                     }
    159                } else {
    160                     goodNetwork = 0; // not a network
    161                }
    162                break;
    163 
    164           default:
    165                goodNetwork = 0; // not a network
    166                break;
    167           }
    168           if (!goodNetwork)
    169                break;
    170      }
    171      if (!goodNetwork) {
    172           delete [] indices_;
    173           // put in message
    174           printf("Not a network - can test if indices_ null\n");
    175           indices_ = NULL;
    176           numberRows_ = 0;
    177           numberColumns_ = 0;
    178      } else {
    179           numberRows_ ++; //  correct
    180           trueNetwork_ = goodNetwork > 0;
    181      }
     69ClpNetworkMatrix::ClpNetworkMatrix(const ClpNetworkMatrix &rhs)
     70  : ClpMatrixBase(rhs)
     71{
     72  matrix_ = NULL;
     73  lengths_ = NULL;
     74  indices_ = NULL;
     75  numberRows_ = rhs.numberRows_;
     76  numberColumns_ = rhs.numberColumns_;
     77  trueNetwork_ = rhs.trueNetwork_;
     78  if (numberColumns_) {
     79    indices_ = new int[2 * numberColumns_];
     80    CoinMemcpyN(rhs.indices_, 2 * numberColumns_, indices_);
     81  }
     82  int numberRows = getNumRows();
     83  if (rhs.rhsOffset_ && numberRows) {
     84    rhsOffset_ = ClpCopyOfArray(rhs.rhsOffset_, numberRows);
     85  } else {
     86    rhsOffset_ = NULL;
     87  }
     88}
     89
     90ClpNetworkMatrix::ClpNetworkMatrix(const CoinPackedMatrix &rhs)
     91  : ClpMatrixBase()
     92{
     93  setType(11);
     94  matrix_ = NULL;
     95  lengths_ = NULL;
     96  indices_ = NULL;
     97  int iColumn;
     98  assert(rhs.isColOrdered());
     99  // get matrix data pointers
     100  const int *row = rhs.getIndices();
     101  const CoinBigIndex *columnStart = rhs.getVectorStarts();
     102  const int *columnLength = rhs.getVectorLengths();
     103  const double *elementByColumn = rhs.getElements();
     104  numberColumns_ = rhs.getNumCols();
     105  int goodNetwork = 1;
     106  numberRows_ = -1;
     107  indices_ = new int[2 * numberColumns_];
     108  CoinBigIndex j = 0;
     109  for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
     110    CoinBigIndex k = columnStart[iColumn];
     111    int iRow;
     112    switch (columnLength[iColumn]) {
     113    case 0:
     114      goodNetwork = -1; // not classic network
     115      indices_[j] = -1;
     116      indices_[j + 1] = -1;
     117      break;
     118
     119    case 1:
     120      goodNetwork = -1; // not classic network
     121      if (fabs(elementByColumn[k] - 1.0) < 1.0e-10) {
     122        indices_[j] = -1;
     123        iRow = row[k];
     124        numberRows_ = CoinMax(numberRows_, iRow);
     125        indices_[j + 1] = iRow;
     126      } else if (fabs(elementByColumn[k] + 1.0) < 1.0e-10) {
     127        indices_[j + 1] = -1;
     128        iRow = row[k];
     129        numberRows_ = CoinMax(numberRows_, iRow);
     130        indices_[j] = iRow;
     131      } else {
     132        goodNetwork = 0; // not a network
     133      }
     134      break;
     135
     136    case 2:
     137      if (fabs(elementByColumn[k] - 1.0) < 1.0e-10) {
     138        if (fabs(elementByColumn[k + 1] + 1.0) < 1.0e-10) {
     139          iRow = row[k];
     140          numberRows_ = CoinMax(numberRows_, iRow);
     141          indices_[j + 1] = iRow;
     142          iRow = row[k + 1];
     143          numberRows_ = CoinMax(numberRows_, iRow);
     144          indices_[j] = iRow;
     145        } else {
     146          goodNetwork = 0; // not a network
     147        }
     148      } else if (fabs(elementByColumn[k] + 1.0) < 1.0e-10) {
     149        if (fabs(elementByColumn[k + 1] - 1.0) < 1.0e-10) {
     150          iRow = row[k];
     151          numberRows_ = CoinMax(numberRows_, iRow);
     152          indices_[j] = iRow;
     153          iRow = row[k + 1];
     154          numberRows_ = CoinMax(numberRows_, iRow);
     155          indices_[j + 1] = iRow;
     156        } else {
     157          goodNetwork = 0; // not a network
     158        }
     159      } else {
     160        goodNetwork = 0; // not a network
     161      }
     162      break;
     163
     164    default:
     165      goodNetwork = 0; // not a network
     166      break;
     167    }
     168    if (!goodNetwork)
     169      break;
     170  }
     171  if (!goodNetwork) {
     172    delete[] indices_;
     173    // put in message
     174    printf("Not a network - can test if indices_ null\n");
     175    indices_ = NULL;
     176    numberRows_ = 0;
     177    numberColumns_ = 0;
     178  } else {
     179    numberRows_++; //  correct
     180    trueNetwork_ = goodNetwork > 0;
     181  }
    182182}
    183183
     
    185185// Destructor
    186186//-------------------------------------------------------------------
    187 ClpNetworkMatrix::~ClpNetworkMatrix ()
    188 {
    189      delete matrix_;
    190      delete [] lengths_;
    191      delete [] indices_;
     187ClpNetworkMatrix::~ClpNetworkMatrix()
     188{
     189  delete matrix_;
     190  delete[] lengths_;
     191  delete[] indices_;
    192192}
    193193
     
    196196//-------------------------------------------------------------------
    197197ClpNetworkMatrix &
    198 ClpNetworkMatrix::operator=(const ClpNetworkMatrix& rhs)
    199 {
    200      if (this != &rhs) {
    201           ClpMatrixBase::operator=(rhs);
    202           delete matrix_;
    203           delete [] lengths_;
    204           delete [] indices_;
    205           matrix_ = NULL;
    206           lengths_ = NULL;
    207           indices_ = NULL;
    208           numberRows_ = rhs.numberRows_;
    209           numberColumns_ = rhs.numberColumns_;
    210           trueNetwork_ = rhs.trueNetwork_;
    211           if (numberColumns_) {
    212                indices_ = new int [ 2*numberColumns_];
    213                CoinMemcpyN(rhs.indices_, 2 * numberColumns_, indices_);
    214           }
    215      }
    216      return *this;
     198ClpNetworkMatrix::operator=(const ClpNetworkMatrix &rhs)
     199{
     200  if (this != &rhs) {
     201    ClpMatrixBase::operator=(rhs);
     202    delete matrix_;
     203    delete[] lengths_;
     204    delete[] indices_;
     205    matrix_ = NULL;
     206    lengths_ = NULL;
     207    indices_ = NULL;
     208    numberRows_ = rhs.numberRows_;
     209    numberColumns_ = rhs.numberColumns_;
     210    trueNetwork_ = rhs.trueNetwork_;
     211    if (numberColumns_) {
     212      indices_ = new int[2 * numberColumns_];
     213      CoinMemcpyN(rhs.indices_, 2 * numberColumns_, indices_);
     214    }
     215  }
     216  return *this;
    217217}
    218218//-------------------------------------------------------------------
    219219// Clone
    220220//-------------------------------------------------------------------
    221 ClpMatrixBase * ClpNetworkMatrix::clone() const
    222 {
    223      return new ClpNetworkMatrix(*this);
     221ClpMatrixBase *ClpNetworkMatrix::clone() const
     222{
     223  return new ClpNetworkMatrix(*this);
    224224}
    225225
     
    228228ClpNetworkMatrix::reverseOrderedCopy() const
    229229{
    230      // count number in each row
    231      CoinBigIndex * tempP = new CoinBigIndex [numberRows_];
    232      CoinBigIndex * tempN = new CoinBigIndex [numberRows_];
    233      memset(tempP, 0, numberRows_ * sizeof(CoinBigIndex));
    234      memset(tempN, 0, numberRows_ * sizeof(CoinBigIndex));
    235      CoinBigIndex j = 0;
    236      int i;
    237      for (i = 0; i < numberColumns_; i++, j += 2) {
    238           int iRow = indices_[j];
    239           tempN[iRow]++;
    240           iRow = indices_[j+1];
    241           tempP[iRow]++;
    242      }
    243      int * newIndices = new int [2*numberColumns_];
    244      CoinBigIndex * newP = new CoinBigIndex [numberRows_+1];
    245      CoinBigIndex * newN = new CoinBigIndex[numberRows_];
    246      int iRow;
    247      j = 0;
    248      // do starts
    249      for (iRow = 0; iRow < numberRows_; iRow++) {
    250           newP[iRow] = j;
    251           j += tempP[iRow];
    252           tempP[iRow] = newP[iRow];
    253           newN[iRow] = j;
    254           j += tempN[iRow];
    255           tempN[iRow] = newN[iRow];
    256      }
    257      newP[numberRows_] = j;
    258      j = 0;
    259      for (i = 0; i < numberColumns_; i++, j += 2) {
    260           int iRow = indices_[j];
    261           CoinBigIndex put = tempN[iRow];
    262           newIndices[put++] = i;
    263           tempN[iRow] = put;
    264           iRow = indices_[j+1];
    265           put = tempP[iRow];
    266           newIndices[put++] = i;
    267           tempP[iRow] = put;
    268      }
    269      delete [] tempP;
    270      delete [] tempN;
    271      ClpPlusMinusOneMatrix * newCopy = new ClpPlusMinusOneMatrix();
    272      newCopy->passInCopy(numberRows_, numberColumns_,
    273                          false, newIndices, newP, newN);
    274      return newCopy;
     230  // count number in each row
     231  CoinBigIndex *tempP = new CoinBigIndex[numberRows_];
     232  CoinBigIndex *tempN = new CoinBigIndex[numberRows_];
     233  memset(tempP, 0, numberRows_ * sizeof(CoinBigIndex));
     234  memset(tempN, 0, numberRows_ * sizeof(CoinBigIndex));
     235  CoinBigIndex j = 0;
     236  int i;
     237  for (i = 0; i < numberColumns_; i++, j += 2) {
     238    int iRow = indices_[j];
     239    tempN[iRow]++;
     240    iRow = indices_[j + 1];
     241    tempP[iRow]++;
     242  }
     243  int *newIndices = new int[2 * numberColumns_];
     244  CoinBigIndex *newP = new CoinBigIndex[numberRows_ + 1];
     245  CoinBigIndex *newN = new CoinBigIndex[numberRows_];
     246  int iRow;
     247  j = 0;
     248  // do starts
     249  for (iRow = 0; iRow < numberRows_; iRow++) {
     250    newP[iRow] = j;
     251    j += tempP[iRow];
     252    tempP[iRow] = newP[iRow];
     253    newN[iRow] = j;
     254    j += tempN[iRow];
     255    tempN[iRow] = newN[iRow];
     256  }
     257  newP[numberRows_] = j;
     258  j = 0;
     259  for (i = 0; i < numberColumns_; i++, j += 2) {
     260    int iRow = indices_[j];
     261    CoinBigIndex put = tempN[iRow];
     262    newIndices[put++] = i;
     263    tempN[iRow] = put;
     264    iRow = indices_[j + 1];
     265    put = tempP[iRow];
     266    newIndices[put++] = i;
     267    tempP[iRow] = put;
     268  }
     269  delete[] tempP;
     270  delete[] tempN;
     271  ClpPlusMinusOneMatrix *newCopy = new ClpPlusMinusOneMatrix();
     272  newCopy->passInCopy(numberRows_, numberColumns_,
     273    false, newIndices, newP, newN);
     274  return newCopy;
    275275}
    276276//unscaled versions
    277 void
    278 ClpNetworkMatrix::times(double scalar,
    279                         const double * x, double * y) const
    280 {
    281      int iColumn;
    282      CoinBigIndex j = 0;
    283      if (trueNetwork_) {
    284           for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
    285                double value = scalar * x[iColumn];
    286                if (value) {
    287                     int iRowM = indices_[j];
    288                     int iRowP = indices_[j+1];
    289                     y[iRowM] -= value;
    290                     y[iRowP] += value;
    291                }
    292           }
    293      } else {
    294           // skip negative rows
    295           for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
    296                double value = scalar * x[iColumn];
    297                if (value) {
    298                     int iRowM = indices_[j];
    299                     int iRowP = indices_[j+1];
    300                     if (iRowM >= 0)
    301                          y[iRowM] -= value;
    302                     if (iRowP >= 0)
    303                          y[iRowP] += value;
    304                }
    305           }
    306      }
    307 }
    308 void
    309 ClpNetworkMatrix::transposeTimes(double scalar,
    310                                  const double * x, double * y) const
    311 {
    312      int iColumn;
    313      CoinBigIndex j = 0;
    314      if (trueNetwork_) {
    315           for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
    316                double value = y[iColumn];
    317                int iRowM = indices_[j];
    318                int iRowP = indices_[j+1];
    319                value -= scalar * x[iRowM];
    320                value += scalar * x[iRowP];
    321                y[iColumn] = value;
    322           }
    323      } else {
    324           // skip negative rows
    325           for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
    326                double value = y[iColumn];
    327                int iRowM = indices_[j];
    328                int iRowP = indices_[j+1];
    329                if (iRowM >= 0)
    330                     value -= scalar * x[iRowM];
    331                if (iRowP >= 0)
    332                     value += scalar * x[iRowP];
    333                y[iColumn] = value;
    334           }
    335      }
    336 }
    337 void
    338 ClpNetworkMatrix::times(double scalar,
    339                         const double * x, double * y,
    340                         const double * /*rowScale*/,
    341                         const double * /*columnScale*/) const
    342 {
    343      // we know it is not scaled
    344      times(scalar, x, y);
    345 }
    346 void
    347 ClpNetworkMatrix::transposeTimes( double scalar,
    348                                   const double * x, double * y,
    349                                   const double * /*rowScale*/,
    350                                   const double * /*columnScale*/,
    351                                   double * /*spare*/) const
    352 {
    353      // we know it is not scaled
    354      transposeTimes(scalar, x, y);
     277void ClpNetworkMatrix::times(double scalar,
     278  const double *x, double *y) const
     279{
     280  int iColumn;
     281  CoinBigIndex j = 0;
     282  if (trueNetwork_) {
     283    for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
     284      double value = scalar * x[iColumn];
     285      if (value) {
     286        int iRowM = indices_[j];
     287        int iRowP = indices_[j + 1];
     288        y[iRowM] -= value;
     289        y[iRowP] += value;
     290      }
     291    }
     292  } else {
     293    // skip negative rows
     294    for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
     295      double value = scalar * x[iColumn];
     296      if (value) {
     297        int iRowM = indices_[j];
     298        int iRowP = indices_[j + 1];
     299        if (iRowM >= 0)
     300          y[iRowM] -= value;
     301        if (iRowP >= 0)
     302          y[iRowP] += value;
     303      }
     304    }
     305  }
     306}
     307void ClpNetworkMatrix::transposeTimes(double scalar,
     308  const double *x, double *y) const
     309{
     310  int iColumn;
     311  CoinBigIndex j = 0;
     312  if (trueNetwork_) {
     313    for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
     314      double value = y[iColumn];
     315      int iRowM = indices_[j];
     316      int iRowP = indices_[j + 1];
     317      value -= scalar * x[iRowM];
     318      value += scalar * x[iRowP];
     319      y[iColumn] = value;
     320    }
     321  } else {
     322    // skip negative rows
     323    for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
     324      double value = y[iColumn];
     325      int iRowM = indices_[j];
     326      int iRowP = indices_[j + 1];
     327      if (iRowM >= 0)
     328        value -= scalar * x[iRowM];
     329      if (iRowP >= 0)
     330        value += scalar * x[iRowP];
     331      y[iColumn] = value;
     332    }
     333  }
     334}
     335void ClpNetworkMatrix::times(double scalar,
     336  const double *x, double *y,
     337  const double * /*rowScale*/,
     338  const double * /*columnScale*/) const
     339{
     340  // we know it is not scaled
     341  times(scalar, x, y);
     342}
     343void ClpNetworkMatrix::transposeTimes(double scalar,
     344  const double *x, double *y,
     345  const double * /*rowScale*/,
     346  const double * /*columnScale*/,
     347  double * /*spare*/) const
     348{
     349  // we know it is not scaled
     350  transposeTimes(scalar, x, y);
    355351}
    356352/* Return <code>x * A + y</code> in <code>z</code>.
    357353        Squashes small elements and knows about ClpSimplex */
    358 void
    359 ClpNetworkMatrix::transposeTimes(const ClpSimplex * model, double scalar,
    360                                  const CoinIndexedVector * rowArray,
    361                                  CoinIndexedVector * y,
    362                                  CoinIndexedVector * columnArray) const
    363 {
    364      // we know it is not scaled
    365      columnArray->clear();
    366      double * pi = rowArray->denseVector();
    367      int numberNonZero = 0;
    368      int * index = columnArray->getIndices();
    369      double * array = columnArray->denseVector();
    370      int numberInRowArray = rowArray->getNumElements();
    371      // maybe I need one in OsiSimplex
    372      double zeroTolerance = model->zeroTolerance();
    373      int numberRows = model->numberRows();
     354void ClpNetworkMatrix::transposeTimes(const ClpSimplex *model, double scalar,
     355  const CoinIndexedVector *rowArray,
     356  CoinIndexedVector *y,
     357  CoinIndexedVector *columnArray) const
     358{
     359  // we know it is not scaled
     360  columnArray->clear();
     361  double *pi = rowArray->denseVector();
     362  int numberNonZero = 0;
     363  int *index = columnArray->getIndices();
     364  double *array = columnArray->denseVector();
     365  int numberInRowArray = rowArray->getNumElements();
     366  // maybe I need one in OsiSimplex
     367  double zeroTolerance = model->zeroTolerance();
     368  int numberRows = model->numberRows();
    374369#ifndef NO_RTTI
    375      ClpPlusMinusOneMatrix* rowCopy =
    376           dynamic_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
     370  ClpPlusMinusOneMatrix *rowCopy = dynamic_cast< ClpPlusMinusOneMatrix * >(model->rowCopy());
    377371#else
    378      ClpPlusMinusOneMatrix* rowCopy =
    379           static_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
     372  ClpPlusMinusOneMatrix *rowCopy = static_cast< ClpPlusMinusOneMatrix * >(model->rowCopy());
    380373#endif
    381      bool packed = rowArray->packedMode();
    382      double factor = 0.3;
    383      // We may not want to do by row if there may be cache problems
    384      int numberColumns = model->numberColumns();
    385      // It would be nice to find L2 cache size - for moment 512K
    386      // Be slightly optimistic
    387      if (numberColumns * sizeof(double) > 1000000) {
    388           if (numberRows * 10 < numberColumns)
    389                factor = 0.1;
    390           else if (numberRows * 4 < numberColumns)
    391                factor = 0.15;
    392           else if (numberRows * 2 < numberColumns)
    393                factor = 0.2;
    394           //if (model->numberIterations()%50==0)
    395           //printf("%d nonzero\n",numberInRowArray);
    396      }
    397      if (numberInRowArray > factor * numberRows || !rowCopy) {
    398           // do by column
    399           int iColumn;
    400           assert (!y->getNumElements());
    401           CoinBigIndex j = 0;
    402           if (packed) {
    403                // need to expand pi into y
    404                assert(y->capacity() >= numberRows);
    405                double * piOld = pi;
    406                pi = y->denseVector();
    407                const int * whichRow = rowArray->getIndices();
    408                int i;
    409                // modify pi so can collapse to one loop
    410                for (i = 0; i < numberInRowArray; i++) {
    411                     int iRow = whichRow[i];
    412                     pi[iRow] = scalar * piOld[i];
    413                }
    414                if (trueNetwork_) {
    415                     for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
    416                          double value = 0.0;
    417                          int iRowM = indices_[j];
    418                          int iRowP = indices_[j+1];
    419                          value -= pi[iRowM];
    420                          value += pi[iRowP];
    421                          if (fabs(value) > zeroTolerance) {
    422                               array[numberNonZero] = value;
    423                               index[numberNonZero++] = iColumn;
    424                          }
    425                     }
    426                } else {
    427                     // skip negative rows
    428                     for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
    429                          double value = 0.0;
    430                          int iRowM = indices_[j];
    431                          int iRowP = indices_[j+1];
    432                          if (iRowM >= 0)
    433                               value -= pi[iRowM];
    434                          if (iRowP >= 0)
    435                               value += pi[iRowP];
    436                          if (fabs(value) > zeroTolerance) {
    437                               array[numberNonZero] = value;
    438                               index[numberNonZero++] = iColumn;
    439                          }
    440                     }
    441                }
    442                for (i = 0; i < numberInRowArray; i++) {
    443                     int iRow = whichRow[i];
    444                     pi[iRow] = 0.0;
    445                }
    446           } else {
    447                if (trueNetwork_) {
    448                     for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
    449                          double value = 0.0;
    450                          int iRowM = indices_[j];
    451                          int iRowP = indices_[j+1];
    452                          value -= scalar * pi[iRowM];
    453                          value += scalar * pi[iRowP];
    454                          if (fabs(value) > zeroTolerance) {
    455                               index[numberNonZero++] = iColumn;
    456                               array[iColumn] = value;
    457                          }
    458                     }
    459                } else {
    460                     // skip negative rows
    461                     for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
    462                          double value = 0.0;
    463                          int iRowM = indices_[j];
    464                          int iRowP = indices_[j+1];
    465                          if (iRowM >= 0)
    466                               value -= scalar * pi[iRowM];
    467                          if (iRowP >= 0)
    468                               value += scalar * pi[iRowP];
    469                          if (fabs(value) > zeroTolerance) {
    470                               index[numberNonZero++] = iColumn;
    471                               array[iColumn] = value;
    472                          }
    473                     }
    474                }
     374  bool packed = rowArray->packedMode();
     375  double factor = 0.3;
     376  // We may not want to do by row if there may be cache problems
     377  int numberColumns = model->numberColumns();
     378  // It would be nice to find L2 cache size - for moment 512K
     379  // Be slightly optimistic
     380  if (numberColumns * sizeof(double) > 1000000) {
     381    if (numberRows * 10 < numberColumns)
     382      factor = 0.1;
     383    else if (numberRows * 4 < numberColumns)
     384      factor = 0.15;
     385    else if (numberRows * 2 < numberColumns)
     386      factor = 0.2;
     387    //if (model->numberIterations()%50==0)
     388    //printf("%d nonzero\n",numberInRowArray);
     389  }
     390  if (numberInRowArray > factor * numberRows || !rowCopy) {
     391    // do by column
     392    int iColumn;
     393    assert(!y->getNumElements());
     394    CoinBigIndex j = 0;
     395    if (packed) {
     396      // need to expand pi into y
     397      assert(y->capacity() >= numberRows);
     398      double *piOld = pi;
     399      pi = y->denseVector();
     400      const int *whichRow = rowArray->getIndices();
     401      int i;
     402      // modify pi so can collapse to one loop
     403      for (i = 0; i < numberInRowArray; i++) {
     404        int iRow = whichRow[i];
     405        pi[iRow] = scalar * piOld[i];
     406      }
     407      if (trueNetwork_) {
     408        for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
     409          double value = 0.0;
     410          int iRowM = indices_[j];
     411          int iRowP = indices_[j + 1];
     412          value -= pi[iRowM];
     413          value += pi[iRowP];
     414          if (fabs(value) > zeroTolerance) {
     415            array[numberNonZero] = value;
     416            index[numberNonZero++] = iColumn;
    475417          }
    476           columnArray->setNumElements(numberNonZero);
    477      } else {
    478           // do by row
    479           rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray);
    480      }
     418        }
     419      } else {
     420        // skip negative rows
     421        for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
     422          double value = 0.0;
     423          int iRowM = indices_[j];
     424          int iRowP = indices_[j + 1];
     425          if (iRowM >= 0)
     426            value -= pi[iRowM];
     427          if (iRowP >= 0)
     428            value += pi[iRowP];
     429          if (fabs(value) > zeroTolerance) {
     430            array[numberNonZero] = value;
     431            index[numberNonZero++] = iColumn;
     432          }
     433        }
     434      }
     435      for (i = 0; i < numberInRowArray; i++) {
     436        int iRow = whichRow[i];
     437        pi[iRow] = 0.0;
     438      }
     439    } else {
     440      if (trueNetwork_) {
     441        for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
     442          double value = 0.0;
     443          int iRowM = indices_[j];
     444          int iRowP = indices_[j + 1];
     445          value -= scalar * pi[iRowM];
     446          value += scalar * pi[iRowP];
     447          if (fabs(value) > zeroTolerance) {
     448            index[numberNonZero++] = iColumn;
     449            array[iColumn] = value;
     450          }
     451        }
     452      } else {
     453        // skip negative rows
     454        for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
     455          double value = 0.0;
     456          int iRowM = indices_[j];
     457          int iRowP = indices_[j + 1];
     458          if (iRowM >= 0)
     459            value -= scalar * pi[iRowM];
     460          if (iRowP >= 0)
     461            value += scalar * pi[iRowP];
     462          if (fabs(value) > zeroTolerance) {
     463            index[numberNonZero++] = iColumn;
     464            array[iColumn] = value;
     465          }
     466        }
     467      }
     468    }
     469    columnArray->setNumElements(numberNonZero);
     470  } else {
     471    // do by row
     472    rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray);
     473  }
    481474}
    482475/* Return <code>x *A in <code>z</code> but
    483476   just for indices in y. */
    484 void
    485 ClpNetworkMatrix::subsetTransposeTimes(const ClpSimplex * /*model*/,
    486                                        const CoinIndexedVector * rowArray,
    487                                        const CoinIndexedVector * y,
    488                                        CoinIndexedVector * columnArray) const
    489 {
    490      columnArray->clear();
    491      double * pi = rowArray->denseVector();
    492      double * array = columnArray->denseVector();
    493      int jColumn;
    494      int numberToDo = y->getNumElements();
    495      const int * which = y->getIndices();
    496      assert (!rowArray->packedMode());
    497      columnArray->setPacked();
    498      if (trueNetwork_) {
    499           for (jColumn = 0; jColumn < numberToDo; jColumn++) {
    500                int iColumn = which[jColumn];
    501                double value = 0.0;
    502                CoinBigIndex j = iColumn << 1;
    503                int iRowM = indices_[j];
    504                int iRowP = indices_[j+1];
    505                value -= pi[iRowM];
    506                value += pi[iRowP];
    507                array[jColumn] = value;
    508           }
    509      } else {
    510           // skip negative rows
    511           for (jColumn = 0; jColumn < numberToDo; jColumn++) {
    512                int iColumn = which[jColumn];
    513                double value = 0.0;
    514                CoinBigIndex j = iColumn << 1;
    515                int iRowM = indices_[j];
    516                int iRowP = indices_[j+1];
    517                if (iRowM >= 0)
    518                     value -= pi[iRowM];
    519                if (iRowP >= 0)
    520                     value += pi[iRowP];
    521                array[jColumn] = value;
    522           }
    523      }
     477void ClpNetworkMatrix::subsetTransposeTimes(const ClpSimplex * /*model*/,
     478  const CoinIndexedVector *rowArray,
     479  const CoinIndexedVector *y,
     480  CoinIndexedVector *columnArray) const
     481{
     482  columnArray->clear();
     483  double *pi = rowArray->denseVector();
     484  double *array = columnArray->denseVector();
     485  int jColumn;
     486  int numberToDo = y->getNumElements();
     487  const int *which = y->getIndices();
     488  assert(!rowArray->packedMode());
     489  columnArray->setPacked();
     490  if (trueNetwork_) {
     491    for (jColumn = 0; jColumn < numberToDo; jColumn++) {
     492      int iColumn = which[jColumn];
     493      double value = 0.0;
     494      CoinBigIndex j = iColumn << 1;
     495      int iRowM = indices_[j];
     496      int iRowP = indices_[j + 1];
     497      value -= pi[iRowM];
     498      value += pi[iRowP];
     499      array[jColumn] = value;
     500    }
     501  } else {
     502    // skip negative rows
     503    for (jColumn = 0; jColumn < numberToDo; jColumn++) {
     504      int iColumn = which[jColumn];
     505      double value = 0.0;
     506      CoinBigIndex j = iColumn << 1;
     507      int iRowM = indices_[j];
     508      int iRowP = indices_[j + 1];
     509      if (iRowM >= 0)
     510        value -= pi[iRowM];
     511      if (iRowP >= 0)
     512        value += pi[iRowP];
     513      array[jColumn] = value;
     514    }
     515  }
    524516}
    525517/// returns number of elements in column part of basis,
    526 int
    527 ClpNetworkMatrix::countBasis( const int * whichColumn,
    528                               int & numberColumnBasic)
    529 {
    530      int i;
    531      CoinBigIndex numberElements = 0;
    532      if (trueNetwork_) {
    533           numberElements = 2 * numberColumnBasic;
    534      } else {
    535           for (i = 0; i < numberColumnBasic; i++) {
    536                int iColumn = whichColumn[i];
    537                CoinBigIndex j = iColumn << 1;
    538                int iRowM = indices_[j];
    539                int iRowP = indices_[j+1];
    540                if (iRowM >= 0)
    541                     numberElements ++;
    542                if (iRowP >= 0)
    543                     numberElements ++;
    544           }
    545      }
     518int ClpNetworkMatrix::countBasis(const int *whichColumn,
     519  int &numberColumnBasic)
     520{
     521  int i;
     522  CoinBigIndex numberElements = 0;
     523  if (trueNetwork_) {
     524    numberElements = 2 * numberColumnBasic;
     525  } else {
     526    for (i = 0; i < numberColumnBasic; i++) {
     527      int iColumn = whichColumn[i];
     528      CoinBigIndex j = iColumn << 1;
     529      int iRowM = indices_[j];
     530      int iRowP = indices_[j + 1];
     531      if (iRowM >= 0)
     532        numberElements++;
     533      if (iRowP >= 0)
     534        numberElements++;
     535    }
     536  }
    546537#if COIN_BIG_INDEX
    547      if (numberElements>COIN_INT_MAX) {
    548        printf("Factorization too large\n");
    549        abort();
    550      }
     538  if (numberElements > COIN_INT_MAX) {
     539    printf("Factorization too large\n");
     540    abort();
     541  }
    551542#endif
    552      return static_cast<int>(numberElements);
    553 }
    554 void
    555 ClpNetworkMatrix::fillBasis(ClpSimplex * /*model*/,
    556                             const int * whichColumn,
    557                             int & numberColumnBasic,
    558                             int * indexRowU, int * start,
    559                             int * rowCount, int * columnCount,
    560                             CoinFactorizationDouble * elementU)
    561 {
    562      int i;
    563      CoinBigIndex numberElements = start[0];
    564      if (trueNetwork_) {
    565           for (i = 0; i < numberColumnBasic; i++) {
    566                int iColumn = whichColumn[i];
    567                CoinBigIndex j = iColumn << 1;
    568                int iRowM = indices_[j];
    569                int iRowP = indices_[j+1];
    570                indexRowU[numberElements] = iRowM;
    571                rowCount[iRowM]++;
    572                elementU[numberElements] = -1.0;
    573                indexRowU[numberElements+1] = iRowP;
    574                rowCount[iRowP]++;
    575                elementU[numberElements+1] = 1.0;
    576                numberElements += 2;
    577                start[i+1] = static_cast<int>(numberElements);
    578                columnCount[i] = 2;
    579           }
    580      } else {
    581           for (i = 0; i < numberColumnBasic; i++) {
    582                int iColumn = whichColumn[i];
    583                CoinBigIndex j = iColumn << 1;
    584                int iRowM = indices_[j];
    585                int iRowP = indices_[j+1];
    586                if (iRowM >= 0) {
    587                     indexRowU[numberElements] = iRowM;
    588                     rowCount[iRowM]++;
    589                     elementU[numberElements++] = -1.0;
    590                }
    591                if (iRowP >= 0) {
    592                     indexRowU[numberElements] = iRowP;
    593                     rowCount[iRowP]++;
    594                     elementU[numberElements++] = 1.0;
    595                }
    596                start[i+1] = static_cast<int>(numberElements);
    597                columnCount[i] = static_cast<int>(numberElements) - start[i];
    598           }
    599      }
     543  return static_cast< int >(numberElements);
     544}
     545void ClpNetworkMatrix::fillBasis(ClpSimplex * /*model*/,
     546  const int *whichColumn,
     547  int &numberColumnBasic,
     548  int *indexRowU, int *start,
     549  int *rowCount, int *columnCount,
     550  CoinFactorizationDouble *elementU)
     551{
     552  int i;
     553  CoinBigIndex numberElements = start[0];
     554  if (trueNetwork_) {
     555    for (i = 0; i < numberColumnBasic; i++) {
     556      int iColumn = whichColumn[i];
     557      CoinBigIndex j = iColumn << 1;
     558      int iRowM = indices_[j];
     559      int iRowP = indices_[j + 1];
     560      indexRowU[numberElements] = iRowM;
     561      rowCount[iRowM]++;
     562      elementU[numberElements] = -1.0;
     563      indexRowU[numberElements + 1] = iRowP;
     564      rowCount[iRowP]++;
     565      elementU[numberElements + 1] = 1.0;
     566      numberElements += 2;
     567      start[i + 1] = static_cast< int >(numberElements);
     568      columnCount[i] = 2;
     569    }
     570  } else {
     571    for (i = 0; i < numberColumnBasic; i++) {
     572      int iColumn = whichColumn[i];
     573      CoinBigIndex j = iColumn << 1;
     574      int iRowM = indices_[j];
     575      int iRowP = indices_[j + 1];
     576      if (iRowM >= 0) {
     577        indexRowU[numberElements] = iRowM;
     578        rowCount[iRowM]++;
     579        elementU[numberElements++] = -1.0;
     580      }
     581      if (iRowP >= 0) {
     582        indexRowU[numberElements] = iRowP;
     583        rowCount[iRowP]++;
     584        elementU[numberElements++] = 1.0;
     585      }
     586      start[i + 1] = static_cast< int >(numberElements);
     587      columnCount[i] = static_cast< int >(numberElements) - start[i];
     588    }
     589  }
    600590#if COIN_BIG_INDEX
    601      if (numberElements>COIN_INT_MAX)
    602        abort();
     591  if (numberElements > COIN_INT_MAX)
     592    abort();
    603593#endif
    604594}
    605595/* Unpacks a column into an CoinIndexedvector
    606596 */
    607 void
    608 ClpNetworkMatrix::unpack(const ClpSimplex * /*model*/, CoinIndexedVector * rowArray,
    609                          int iColumn) const
    610 {
    611      CoinBigIndex j = iColumn << 1;
    612      int iRowM = indices_[j];
    613      int iRowP = indices_[j+1];
    614      if (iRowM >= 0)
    615           rowArray->add(iRowM, -1.0);
    616      if (iRowP >= 0)
    617           rowArray->add(iRowP, 1.0);
     597void ClpNetworkMatrix::unpack(const ClpSimplex * /*model*/, CoinIndexedVector *rowArray,
     598  int iColumn) const
     599{
     600  CoinBigIndex j = iColumn << 1;
     601  int iRowM = indices_[j];
     602  int iRowP = indices_[j + 1];
     603  if (iRowM >= 0)
     604    rowArray->add(iRowM, -1.0);
     605  if (iRowP >= 0)
     606    rowArray->add(iRowP, 1.0);
    618607}
    619608/* Unpacks a column into an CoinIndexedvector
     
    621610Note that model is NOT const.  Bounds and objective could
    622611be modified if doing column generation (just for this variable) */
    623 void
    624 ClpNetworkMatrix::unpackPacked(ClpSimplex * /*model*/,
    625                                CoinIndexedVector * rowArray,
    626                                int iColumn) const
    627 {
    628      int * index = rowArray->getIndices();
    629      double * array = rowArray->denseVector();
    630      int number = 0;
    631      CoinBigIndex j = iColumn << 1;
    632      int iRowM = indices_[j];
    633      int iRowP = indices_[j+1];
    634      if (iRowM >= 0) {
    635           array[number] = -1.0;
    636           index[number++] = iRowM;
    637      }
    638      if (iRowP >= 0) {
    639           array[number] = 1.0;
    640           index[number++] = iRowP;
    641      }
    642      rowArray->setNumElements(number);
    643      rowArray->setPackedMode(true);
     612void ClpNetworkMatrix::unpackPacked(ClpSimplex * /*model*/,
     613  CoinIndexedVector *rowArray,
     614  int iColumn) const
     615{
     616  int *index = rowArray->getIndices();
     617  double *array = rowArray->denseVector();
     618  int number = 0;
     619  CoinBigIndex j = iColumn << 1;
     620  int iRowM = indices_[j];
     621  int iRowP = indices_[j + 1];
     622  if (iRowM >= 0) {
     623    array[number] = -1.0;
     624    index[number++] = iRowM;
     625  }
     626  if (iRowP >= 0) {
     627    array[number] = 1.0;
     628    index[number++] = iRowP;
     629  }
     630  rowArray->setNumElements(number);
     631  rowArray->setPackedMode(true);
    644632}
    645633/* Adds multiple of a column into an CoinIndexedvector
    646634      You can use quickAdd to add to vector */
    647 void
    648 ClpNetworkMatrix::add(const ClpSimplex * /*model*/, CoinIndexedVector * rowArray,
    649                       int iColumn, double multiplier) const
    650 {
    651      CoinBigIndex j = iColumn << 1;
    652      int iRowM = indices_[j];
    653      int iRowP = indices_[j+1];
    654      if (iRowM >= 0)
    655           rowArray->quickAdd(iRowM, -multiplier);
    656      if (iRowP >= 0)
    657           rowArray->quickAdd(iRowP, multiplier);
     635void ClpNetworkMatrix::add(const ClpSimplex * /*model*/, CoinIndexedVector *rowArray,
     636  int iColumn, double multiplier) const
     637{
     638  CoinBigIndex j = iColumn << 1;
     639  int iRowM = indices_[j];
     640  int iRowP = indices_[j + 1];
     641  if (iRowM >= 0)
     642    rowArray->quickAdd(iRowM, -multiplier);
     643  if (iRowP >= 0)
     644    rowArray->quickAdd(iRowP, multiplier);
    658645}
    659646/* Adds multiple of a column into an array */
    660 void
    661 ClpNetworkMatrix::add(const ClpSimplex * /*model*/, double * array,
    662                       int iColumn, double multiplier) const
    663 {
    664      CoinBigIndex j = iColumn << 1;
    665      int iRowM = indices_[j];
    666      int iRowP = indices_[j+1];
    667      if (iRowM >= 0)
    668           array[iRowM] -= multiplier;
    669      if (iRowP >= 0)
    670           array[iRowP] += multiplier;
     647void ClpNetworkMatrix::add(const ClpSimplex * /*model*/, double *array,
     648  int iColumn, double multiplier) const
     649{
     650  CoinBigIndex j = iColumn << 1;
     651  int iRowM = indices_[j];
     652  int iRowP = indices_[j + 1];
     653  if (iRowM >= 0)
     654    array[iRowM] -= multiplier;
     655  if (iRowP >= 0)
     656    array[iRowP] += multiplier;
    671657}
    672658
     
    675661ClpNetworkMatrix::getPackedMatrix() const
    676662{
    677      if (!matrix_) {
    678           assert (trueNetwork_); // fix later
    679           int numberElements = 2 * numberColumns_;
    680           double * elements = new double [numberElements];
    681           CoinBigIndex i;
    682           for (i = 0; i < 2 * numberColumns_; i += 2) {
    683                elements[i] = -1.0;
    684                elements[i+1] = 1.0;
    685           }
    686           CoinBigIndex * starts = new CoinBigIndex [numberColumns_+1];
    687           for (i = 0; i < numberColumns_ + 1; i++) {
    688                starts[i] = 2 * i;
    689           }
    690           // use assignMatrix to save space
    691           delete [] lengths_;
    692           lengths_ = NULL;
    693           matrix_ = new CoinPackedMatrix();
    694           int * indices = CoinCopyOfArray(indices_, 2 * numberColumns_);
    695           matrix_->assignMatrix(true, numberRows_, numberColumns_,
    696                                 getNumElements(),
    697                                 elements, indices,
    698                                 starts, lengths_);
    699           assert(!elements);
    700           assert(!starts);
    701           assert (!indices);
    702           assert (!lengths_);
    703      }
    704      return matrix_;
     663  if (!matrix_) {
     664    assert(trueNetwork_); // fix later
     665    int numberElements = 2 * numberColumns_;
     666    double *elements = new double[numberElements];
     667    CoinBigIndex i;
     668    for (i = 0; i < 2 * numberColumns_; i += 2) {
     669      elements[i] = -1.0;
     670      elements[i + 1] = 1.0;
     671    }
     672    CoinBigIndex *starts = new CoinBigIndex[numberColumns_ + 1];
     673    for (i = 0; i < numberColumns_ + 1; i++) {
     674      starts[i] = 2 * i;
     675    }
     676    // use assignMatrix to save space
     677    delete[] lengths_;
     678    lengths_ = NULL;
     679    matrix_ = new CoinPackedMatrix();
     680    int *indices = CoinCopyOfArray(indices_, 2 * numberColumns_);
     681    matrix_->assignMatrix(true, numberRows_, numberColumns_,
     682      getNumElements(),
     683      elements, indices,
     684      starts, lengths_);
     685    assert(!elements);
     686    assert(!starts);
     687    assert(!indices);
     688    assert(!lengths_);
     689  }
     690  return matrix_;
    705691}
    706692/* A vector containing the elements in the packed matrix. Note that there
     
    711697ClpNetworkMatrix::getElements() const
    712698{
    713      if (!matrix_)
    714           getPackedMatrix();
    715      return matrix_->getElements();
     699  if (!matrix_)
     700    getPackedMatrix();
     701  return matrix_->getElements();
    716702}
    717703
     
    719705ClpNetworkMatrix::getVectorStarts() const
    720706{
    721      if (!matrix_)
    722           getPackedMatrix();
    723      return matrix_->getVectorStarts();
     707  if (!matrix_)
     708    getPackedMatrix();
     709  return matrix_->getVectorStarts();
    724710}
    725711/* The lengths of the major-dimension vectors. */
     
    727713ClpNetworkMatrix::getVectorLengths() const
    728714{
    729      assert (trueNetwork_); // fix later
    730      if (!lengths_) {
    731           lengths_ = new int [numberColumns_];
    732           int i;
    733           for (i = 0; i < numberColumns_; i++) {
    734                lengths_[i] = 2;
    735           }
    736      }
    737      return lengths_;
     715  assert(trueNetwork_); // fix later
     716  if (!lengths_) {
     717    lengths_ = new int[numberColumns_];
     718    int i;
     719    for (i = 0; i < numberColumns_; i++) {
     720      lengths_[i] = 2;
     721    }
     722  }
     723  return lengths_;
    738724}
    739725/* Delete the columns whose indices are listed in <code>indDel</code>. */
    740 void ClpNetworkMatrix::deleteCols(const int numDel, const int * indDel)
    741 {
    742      assert (trueNetwork_);
    743      int iColumn;
    744      int numberBad = 0;
    745      // Use array to make sure we can have duplicates
    746      char * which = new char[numberColumns_];
    747      memset(which, 0, numberColumns_);
    748      int nDuplicate = 0;
    749      for (iColumn = 0; iColumn < numDel; iColumn++) {
    750           int jColumn = indDel[iColumn];
    751           if (jColumn < 0 || jColumn >= numberColumns_) {
    752                numberBad++;
    753           } else {
    754                if (which[jColumn])
    755                     nDuplicate++;
    756                else
    757                     which[jColumn] = 1;
    758           }
    759      }
    760      if (numberBad)
    761           throw CoinError("Indices out of range", "deleteCols", "ClpNetworkMatrix");
    762      int newNumber = numberColumns_ - numDel + nDuplicate;
    763      // Get rid of temporary arrays
    764      delete [] lengths_;
    765      lengths_ = NULL;
    766      delete matrix_;
    767      matrix_ = NULL;
    768      int newSize = 2 * newNumber;
    769      int * newIndices = new int [newSize];
    770      newSize = 0;
    771      for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    772           if (!which[iColumn]) {
    773                CoinBigIndex start;
    774                CoinBigIndex i;
    775                start = 2 * iColumn;
    776                for (i = start; i < start + 2; i++)
    777                     newIndices[newSize++] = indices_[i];
    778           }
    779      }
    780      delete [] which;
    781      delete [] indices_;
    782      indices_ = newIndices;
    783      numberColumns_ = newNumber;
     726void ClpNetworkMatrix::deleteCols(const int numDel, const int *indDel)
     727{
     728  assert(trueNetwork_);
     729  int iColumn;
     730  int numberBad = 0;
     731  // Use array to make sure we can have duplicates
     732  char *which = new char[numberColumns_];
     733  memset(which, 0, numberColumns_);
     734  int nDuplicate = 0;
     735  for (iColumn = 0; iColumn < numDel; iColumn++) {
     736    int jColumn = indDel[iColumn];
     737    if (jColumn < 0 || jColumn >= numberColumns_) {
     738      numberBad++;
     739    } else {
     740      if (which[jColumn])
     741        nDuplicate++;
     742      else
     743        which[jColumn] = 1;
     744    }
     745  }
     746  if (numberBad)
     747    throw CoinError("Indices out of range", "deleteCols", "ClpNetworkMatrix");
     748  int newNumber = numberColumns_ - numDel + nDuplicate;
     749  // Get rid of temporary arrays
     750  delete[] lengths_;
     751  lengths_ = NULL;
     752  delete matrix_;
     753  matrix_ = NULL;
     754  int newSize = 2 * newNumber;
     755  int *newIndices = new int[newSize];
     756  newSize = 0;
     757  for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     758    if (!which[iColumn]) {
     759      CoinBigIndex start;
     760      CoinBigIndex i;
     761      start = 2 * iColumn;
     762      for (i = start; i < start + 2; i++)
     763        newIndices[newSize++] = indices_[i];
     764    }
     765  }
     766  delete[] which;
     767  delete[] indices_;
     768  indices_ = newIndices;
     769  numberColumns_ = newNumber;
    784770}
    785771/* Delete the rows whose indices are listed in <code>indDel</code>. */
    786 void ClpNetworkMatrix::deleteRows(const int numDel, const int * indDel)
    787 {
    788      int iRow;
    789      int numberBad = 0;
    790      // Use array to make sure we can have duplicates
    791      int * which = new int [numberRows_];
    792      memset(which, 0, numberRows_ * sizeof(int));
    793      for (iRow = 0; iRow < numDel; iRow++) {
    794           int jRow = indDel[iRow];
    795           if (jRow < 0 || jRow >= numberRows_) {
    796                numberBad++;
    797           } else {
    798                which[jRow] = 1;
    799           }
    800      }
    801      if (numberBad)
    802           throw CoinError("Indices out of range", "deleteRows", "ClpNetworkMatrix");
    803      // Only valid of all columns have 0 entries
    804      int iColumn;
    805      for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    806           CoinBigIndex start;
    807           CoinBigIndex i;
    808           start = 2 * iColumn;
    809           for (i = start; i < start + 2; i++) {
    810                int iRow = indices_[i];
    811                if (which[iRow])
    812                     numberBad++;
    813           }
    814      }
    815      if (numberBad)
    816           throw CoinError("Row has entries", "deleteRows", "ClpNetworkMatrix");
    817      int newNumber = 0;
    818      for (iRow = 0; iRow < numberRows_; iRow++) {
    819           if (!which[iRow])
    820                which[iRow] = newNumber++;
    821           else
    822                which[iRow] = -1;
    823      }
    824      for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    825           CoinBigIndex start;
    826           CoinBigIndex i;
    827           start = 2 * iColumn;
    828           for (i = start; i < start + 2; i++) {
    829                int iRow = indices_[i];
    830                indices_[i] = which[iRow];
    831           }
    832      }
    833      delete [] which;
    834      numberRows_ = newNumber;
     772void ClpNetworkMatrix::deleteRows(const int numDel, const int *indDel)
     773{
     774  int iRow;
     775  int numberBad = 0;
     776  // Use array to make sure we can have duplicates
     777  int *which = new int[numberRows_];
     778  memset(which, 0, numberRows_ * sizeof(int));
     779  for (iRow = 0; iRow < numDel; iRow++) {
     780    int jRow = indDel[iRow];
     781    if (jRow < 0 || jRow >= numberRows_) {
     782      numberBad++;
     783    } else {
     784      which[jRow] = 1;
     785    }
     786  }
     787  if (numberBad)
     788    throw CoinError("Indices out of range", "deleteRows", "ClpNetworkMatrix");
     789  // Only valid of all columns have 0 entries
     790  int iColumn;
     791  for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     792    CoinBigIndex start;
     793    CoinBigIndex i;
     794    start = 2 * iColumn;
     795    for (i = start; i < start + 2; i++) {
     796      int iRow = indices_[i];
     797      if (which[iRow])
     798        numberBad++;
     799    }
     800  }
     801  if (numberBad)
     802    throw CoinError("Row has entries", "deleteRows", "ClpNetworkMatrix");
     803  int newNumber = 0;
     804  for (iRow = 0; iRow < numberRows_; iRow++) {
     805    if (!which[iRow])
     806      which[iRow] = newNumber++;
     807    else
     808      which[iRow] = -1;
     809  }
     810  for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     811    CoinBigIndex start;
     812    CoinBigIndex i;
     813    start = 2 * iColumn;
     814    for (i = start; i < start + 2; i++) {
     815      int iRow = indices_[i];
     816      indices_[i] = which[iRow];
     817    }
     818  }
     819  delete[] which;
     820  numberRows_ = newNumber;
    835821}
    836822/* Given positive integer weights for each row fills in sum of weights
     
    839825*/
    840826CoinBigIndex *
    841 ClpNetworkMatrix::dubiousWeights(const ClpSimplex * model, int * inputWeights) const
    842 {
    843      int numberRows = model->numberRows();
    844      int numberColumns = model->numberColumns();
    845      int number = numberRows + numberColumns;
    846      CoinBigIndex * weights = new CoinBigIndex[number];
    847      int i;
    848      for (i = 0; i < numberColumns; i++) {
    849           CoinBigIndex j = i << 1;
    850           CoinBigIndex count = 0;
    851           int iRowM = indices_[j];
    852           int iRowP = indices_[j+1];
    853           if (iRowM >= 0) {
    854                count += inputWeights[iRowM];
    855           }
    856           if (iRowP >= 0) {
    857                count += inputWeights[iRowP];
    858           }
    859           weights[i] = count;
    860      }
    861      for (i = 0; i < numberRows; i++) {
    862           weights[i+numberColumns] = inputWeights[i];
    863      }
    864      return weights;
     827ClpNetworkMatrix::dubiousWeights(const ClpSimplex *model, int *inputWeights) const
     828{
     829  int numberRows = model->numberRows();
     830  int numberColumns = model->numberColumns();
     831  int number = numberRows + numberColumns;
     832  CoinBigIndex *weights = new CoinBigIndex[number];
     833  int i;
     834  for (i = 0; i < numberColumns; i++) {
     835    CoinBigIndex j = i << 1;
     836    CoinBigIndex count = 0;
     837    int iRowM = indices_[j];
     838    int iRowP = indices_[j + 1];
     839    if (iRowM >= 0) {
     840      count += inputWeights[iRowM];
     841    }
     842    if (iRowP >= 0) {
     843      count += inputWeights[iRowP];
     844    }
     845    weights[i] = count;
     846  }
     847  for (i = 0; i < numberRows; i++) {
     848    weights[i + numberColumns] = inputWeights[i];
     849  }
     850  return weights;
    865851}
    866852/* Returns largest and smallest elements of both signs.
    867853   Largest refers to largest absolute value.
    868854*/
    869 void
    870 ClpNetworkMatrix::rangeOfElements(double & smallestNegative, double & largestNegative,
    871                                   double & smallestPositive, double & largestPositive)
    872 {
    873      smallestNegative = -1.0;
    874      largestNegative = -1.0;
    875      smallestPositive = 1.0;
    876      largestPositive = 1.0;
     855void ClpNetworkMatrix::rangeOfElements(double &smallestNegative, double &largestNegative,
     856  double &smallestPositive, double &largestPositive)
     857{
     858  smallestNegative = -1.0;
     859  largestNegative = -1.0;
     860  smallestPositive = 1.0;
     861  largestPositive = 1.0;
    877862}
    878863// Says whether it can do partial pricing
    879 bool
    880 ClpNetworkMatrix::canDoPartialPricing() const
    881 {
    882      return true;
     864bool ClpNetworkMatrix::canDoPartialPricing() const
     865{
     866  return true;
    883867}
    884868// Partial pricing
    885 void
    886 ClpNetworkMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
    887                                  int & bestSequence, int & numberWanted)
    888 {
    889      numberWanted = currentWanted_;
    890      int j;
    891      int start = static_cast<int> (startFraction * numberColumns_);
    892      int end = CoinMin(static_cast<int> (endFraction * numberColumns_ + 1), numberColumns_);
    893      double tolerance = model->currentDualTolerance();
    894      double * reducedCost = model->djRegion();
    895      const double * duals = model->dualRowSolution();
    896      const double * cost = model->costRegion();
    897      double bestDj;
    898      if (bestSequence >= 0)
    899           bestDj = fabs(reducedCost[bestSequence]);
    900      else
    901           bestDj = tolerance;
    902      int sequenceOut = model->sequenceOut();
    903      int saveSequence = bestSequence;
    904      if (!trueNetwork_) {
    905           // Not true network
    906           int iSequence;
    907           for (iSequence = start; iSequence < end; iSequence++) {
    908                if (iSequence != sequenceOut) {
    909                     double value;
    910                     int iRowM, iRowP;
    911                     ClpSimplex::Status status = model->getStatus(iSequence);
    912 
    913                     switch(status) {
    914 
    915                     case ClpSimplex::basic:
    916                     case ClpSimplex::isFixed:
    917                          break;
    918                     case ClpSimplex::isFree:
    919                     case ClpSimplex::superBasic:
    920                          value = cost[iSequence];
    921                          j = iSequence << 1;
    922                          // skip negative rows
    923                          iRowM = indices_[j];
    924                          iRowP = indices_[j+1];
    925                          if (iRowM >= 0)
    926                               value += duals[iRowM];
    927                          if (iRowP >= 0)
    928                               value -= duals[iRowP];
    929                          value = fabs(value);
    930                          if (value > FREE_ACCEPT * tolerance) {
    931                               numberWanted--;
    932                               // we are going to bias towards free (but only if reasonable)
    933                               value *= FREE_BIAS;
    934                               if (value > bestDj) {
    935                                    // check flagged variable and correct dj
    936                                    if (!model->flagged(iSequence)) {
    937                                         bestDj = value;
    938                                         bestSequence = iSequence;
    939                                    } else {
    940                                         // just to make sure we don't exit before got something
    941                                         numberWanted++;
    942                                    }
    943                               }
    944                          }
    945                          break;
    946                     case ClpSimplex::atUpperBound:
    947                          value = cost[iSequence];
    948                          j = iSequence << 1;
    949                          // skip negative rows
    950                          iRowM = indices_[j];
    951                          iRowP = indices_[j+1];
    952                          if (iRowM >= 0)
    953                               value += duals[iRowM];
    954                          if (iRowP >= 0)
    955                               value -= duals[iRowP];
    956                          if (value > tolerance) {
    957                               numberWanted--;
    958                               if (value > bestDj) {
    959                                    // check flagged variable and correct dj
    960                                    if (!model->flagged(iSequence)) {
    961                                         bestDj = value;
    962                                         bestSequence = iSequence;
    963                                    } else {
    964                                         // just to make sure we don't exit before got something
    965                                         numberWanted++;
    966                                    }
    967                               }
    968                          }
    969                          break;
    970                     case ClpSimplex::atLowerBound:
    971                          value = cost[iSequence];
    972                          j = iSequence << 1;
    973                          // skip negative rows
    974                          iRowM = indices_[j];
    975                          iRowP = indices_[j+1];
    976                          if (iRowM >= 0)
    977                               value += duals[iRowM];
    978                          if (iRowP >= 0)
    979                               value -= duals[iRowP];
    980                          value = -value;
    981                          if (value > tolerance) {
    982                               numberWanted--;
    983                               if (value > bestDj) {
    984                                    // check flagged variable and correct dj
    985                                    if (!model->flagged(iSequence)) {
    986                                         bestDj = value;
    987                                         bestSequence = iSequence;
    988                                    } else {
    989                                         // just to make sure we don't exit before got something
    990                                         numberWanted++;
    991                                    }
    992                               }
    993                          }
    994                          break;
    995                     }
    996                }
    997                if (!numberWanted)
    998                     break;
     869void ClpNetworkMatrix::partialPricing(ClpSimplex *model, double startFraction, double endFraction,
     870  int &bestSequence, int &numberWanted)
     871{
     872  numberWanted = currentWanted_;
     873  int j;
     874  int start = static_cast< int >(startFraction * numberColumns_);
     875  int end = CoinMin(static_cast< int >(endFraction * numberColumns_ + 1), numberColumns_);
     876  double tolerance = model->currentDualTolerance();
     877  double *reducedCost = model->djRegion();
     878  const double *duals = model->dualRowSolution();
     879  const double *cost = model->costRegion();
     880  double bestDj;
     881  if (bestSequence >= 0)
     882    bestDj = fabs(reducedCost[bestSequence]);
     883  else
     884    bestDj = tolerance;
     885  int sequenceOut = model->sequenceOut();
     886  int saveSequence = bestSequence;
     887  if (!trueNetwork_) {
     888    // Not true network
     889    int iSequence;
     890    for (iSequence = start; iSequence < end; iSequence++) {
     891      if (iSequence != sequenceOut) {
     892        double value;
     893        int iRowM, iRowP;
     894        ClpSimplex::Status status = model->getStatus(iSequence);
     895
     896        switch (status) {
     897
     898        case ClpSimplex::basic:
     899        case ClpSimplex::isFixed:
     900          break;
     901        case ClpSimplex::isFree:
     902        case ClpSimplex::superBasic:
     903          value = cost[iSequence];
     904          j = iSequence << 1;
     905          // skip negative rows
     906          iRowM = indices_[j];
     907          iRowP = indices_[j + 1];
     908          if (iRowM >= 0)
     909            value += duals[iRowM];
     910          if (iRowP >= 0)
     911            value -= duals[iRowP];
     912          value = fabs(value);
     913          if (value > FREE_ACCEPT * tolerance) {
     914            numberWanted--;
     915            // we are going to bias towards free (but only if reasonable)
     916            value *= FREE_BIAS;
     917            if (value > bestDj) {
     918              // check flagged variable and correct dj
     919              if (!model->flagged(iSequence)) {
     920                bestDj = value;
     921                bestSequence = iSequence;
     922              } else {
     923                // just to make sure we don't exit before got something
     924                numberWanted++;
     925              }
     926            }
    999927          }
    1000           if (bestSequence != saveSequence) {
    1001                // recompute dj
    1002                double value = cost[bestSequence];
    1003                j = bestSequence << 1;
    1004                // skip negative rows
    1005                int iRowM = indices_[j];
    1006                int iRowP = indices_[j+1];
    1007                if (iRowM >= 0)
    1008                     value += duals[iRowM];
    1009                if (iRowP >= 0)
    1010                     value -= duals[iRowP];
    1011                reducedCost[bestSequence] = value;
    1012                savedBestSequence_ = bestSequence;
    1013                savedBestDj_ = reducedCost[savedBestSequence_];
     928          break;
     929        case ClpSimplex::atUpperBound:
     930          value = cost[iSequence];
     931          j = iSequence << 1;
     932          // skip negative rows
     933          iRowM = indices_[j];
     934          iRowP = indices_[j + 1];
     935          if (iRowM >= 0)
     936            value += duals[iRowM];
     937          if (iRowP >= 0)
     938            value -= duals[iRowP];
     939          if (value > tolerance) {
     940            numberWanted--;
     941            if (value > bestDj) {
     942              // check flagged variable and correct dj
     943              if (!model->flagged(iSequence)) {
     944                bestDj = value;
     945                bestSequence = iSequence;
     946              } else {
     947                // just to make sure we don't exit before got something
     948                numberWanted++;
     949              }
     950            }
    1014951          }
    1015      } else {
    1016           // true network
    1017           int iSequence;
    1018           for (iSequence = start; iSequence < end; iSequence++) {
    1019                if (iSequence != sequenceOut) {
    1020                     double value;
    1021                     int iRowM, iRowP;
    1022                     ClpSimplex::Status status = model->getStatus(iSequence);
    1023 
    1024                     switch(status) {
    1025 
    1026                     case ClpSimplex::basic:
    1027                     case ClpSimplex::isFixed:
    1028                          break;
    1029                     case ClpSimplex::isFree:
    1030                     case ClpSimplex::superBasic:
    1031                          value = cost[iSequence];
    1032                          j = iSequence << 1;
    1033                          iRowM = indices_[j];
    1034                          iRowP = indices_[j+1];
    1035                          value += duals[iRowM];
    1036                          value -= duals[iRowP];
    1037                          value = fabs(value);
    1038                          if (value > FREE_ACCEPT * tolerance) {
    1039                               numberWanted--;
    1040                               // we are going to bias towards free (but only if reasonable)
    1041                               value *= FREE_BIAS;
    1042                               if (value > bestDj) {
    1043                                    // check flagged variable and correct dj
    1044                                    if (!model->flagged(iSequence)) {
    1045                                         bestDj = value;
    1046                                         bestSequence = iSequence;
    1047                                    } else {
    1048                                         // just to make sure we don't exit before got something
    1049                                         numberWanted++;
    1050                                    }
    1051                               }
    1052                          }
    1053                          break;
    1054                     case ClpSimplex::atUpperBound:
    1055                          value = cost[iSequence];
    1056                          j = iSequence << 1;
    1057                          iRowM = indices_[j];
    1058                          iRowP = indices_[j+1];
    1059                          value += duals[iRowM];
    1060                          value -= duals[iRowP];
    1061                          if (value > tolerance) {
    1062                               numberWanted--;
    1063                               if (value > bestDj) {
    1064                                    // check flagged variable and correct dj
    1065                                    if (!model->flagged(iSequence)) {
    1066                                         bestDj = value;
    1067                                         bestSequence = iSequence;
    1068                                    } else {
    1069                                         // just to make sure we don't exit before got something
    1070                                         numberWanted++;
    1071                                    }
    1072                               }
    1073                          }
    1074                          break;
    1075                     case ClpSimplex::atLowerBound:
    1076                          value = cost[iSequence];
    1077                          j = iSequence << 1;
    1078                          iRowM = indices_[j];
    1079                          iRowP = indices_[j+1];
    1080                          value += duals[iRowM];
    1081                          value -= duals[iRowP];
    1082                          value = -value;
    1083                          if (value > tolerance) {
    1084                               numberWanted--;
    1085                               if (value > bestDj) {
    1086                                    // check flagged variable and correct dj
    1087                                    if (!model->flagged(iSequence)) {
    1088                                         bestDj = value;
    1089                                         bestSequence = iSequence;
    1090                                    } else {
    1091                                         // just to make sure we don't exit before got something
    1092                                         numberWanted++;
    1093                                    }
    1094                               }
    1095                          }
    1096                          break;
    1097                     }
    1098                }
    1099                if (!numberWanted)
    1100                     break;
     952          break;
     953        case ClpSimplex::atLowerBound:
     954          value = cost[iSequence];
     955          j = iSequence << 1;
     956          // skip negative rows
     957          iRowM = indices_[j];
     958          iRowP = indices_[j + 1];
     959          if (iRowM >= 0)
     960            value += duals[iRowM];
     961          if (iRowP >= 0)
     962            value -= duals[iRowP];
     963          value = -value;
     964          if (value > tolerance) {
     965            numberWanted--;
     966            if (value > bestDj) {
     967              // check flagged variable and correct dj
     968              if (!model->flagged(iSequence)) {
     969                bestDj = value;
     970                bestSequence = iSequence;
     971              } else {
     972                // just to make sure we don't exit before got something
     973                numberWanted++;
     974              }
     975            }
    1101976          }
    1102           if (bestSequence != saveSequence) {
    1103                // recompute dj
    1104                double value = cost[bestSequence];
    1105                j = bestSequence << 1;
    1106                int iRowM = indices_[j];
    1107                int iRowP = indices_[j+1];
    1108                value += duals[iRowM];
    1109                value -= duals[iRowP];
    1110                reducedCost[bestSequence] = value;
    1111                savedBestSequence_ = bestSequence;
    1112                savedBestDj_ = reducedCost[savedBestSequence_];
     977          break;
     978        }
     979      }
     980      if (!numberWanted)
     981        break;
     982    }
     983    if (bestSequence != saveSequence) {
     984      // recompute dj
     985      double value = cost[bestSequence];
     986      j = bestSequence << 1;
     987      // skip negative rows
     988      int iRowM = indices_[j];
     989      int iRowP = indices_[j + 1];
     990      if (iRowM >= 0)
     991        value += duals[iRowM];
     992      if (iRowP >= 0)
     993        value -= duals[iRowP];
     994      reducedCost[bestSequence] = value;
     995      savedBestSequence_ = bestSequence;
     996      savedBestDj_ = reducedCost[savedBestSequence_];
     997    }
     998  } else {
     999    // true network
     1000    int iSequence;
     1001    for (iSequence = start; iSequence < end; iSequence++) {
     1002      if (iSequence != sequenceOut) {
     1003        double value;
     1004        int iRowM, iRowP;
     1005        ClpSimplex::Status status = model->getStatus(iSequence);
     1006
     1007        switch (status) {
     1008
     1009        case ClpSimplex::basic:
     1010        case ClpSimplex::isFixed:
     1011          break;
     1012        case ClpSimplex::isFree:
     1013        case ClpSimplex::superBasic:
     1014          value = cost[iSequence];
     1015          j = iSequence << 1;
     1016          iRowM = indices_[j];
     1017          iRowP = indices_[j + 1];
     1018          value += duals[iRowM];
     1019          value -= duals[iRowP];
     1020          value = fabs(value);
     1021          if (value > FREE_ACCEPT * tolerance) {
     1022            numberWanted--;
     1023            // we are going to bias towards free (but only if reasonable)
     1024            value *= FREE_BIAS;
     1025            if (value > bestDj) {
     1026              // check flagged variable and correct dj
     1027              if (!model->flagged(iSequence)) {
     1028                bestDj = value;
     1029                bestSequence = iSequence;
     1030              } else {
     1031                // just to make sure we don't exit before got something
     1032                numberWanted++;
     1033              }
     1034            }
    11131035          }
    1114      }
    1115      currentWanted_ = numberWanted;
     1036          break;
     1037        case ClpSimplex::atUpperBound:
     1038          value = cost[iSequence];
     1039          j = iSequence << 1;
     1040          iRowM = indices_[j];
     1041          iRowP = indices_[j + 1];
     1042          value += duals[iRowM];
     1043          value -= duals[iRowP];
     1044          if (value > tolerance) {
     1045            numberWanted--;
     1046            if (value > bestDj) {
     1047              // check flagged variable and correct dj
     1048              if (!model->flagged(iSequence)) {
     1049                bestDj = value;
     1050                bestSequence = iSequence;
     1051              } else {
     1052                // just to make sure we don't exit before got something
     1053                numberWanted++;
     1054              }
     1055            }
     1056          }
     1057          break;
     1058        case ClpSimplex::atLowerBound:
     1059          value = cost[iSequence];
     1060          j = iSequence << 1;
     1061          iRowM = indices_[j];
     1062          iRowP = indices_[j + 1];
     1063          value += duals[iRowM];
     1064          value -= duals[iRowP];
     1065          value = -value;
     1066          if (value > tolerance) {
     1067            numberWanted--;
     1068            if (value > bestDj) {
     1069              // check flagged variable and correct dj
     1070              if (!model->flagged(iSequence)) {
     1071                bestDj = value;
     1072                bestSequence = iSequence;
     1073              } else {
     1074                // just to make sure we don't exit before got something
     1075                numberWanted++;
     1076              }
     1077            }
     1078          }
     1079          break;
     1080        }
     1081      }
     1082      if (!numberWanted)
     1083        break;
     1084    }
     1085    if (bestSequence != saveSequence) {
     1086      // recompute dj
     1087      double value = cost[bestSequence];
     1088      j = bestSequence << 1;
     1089      int iRowM = indices_[j];
     1090      int iRowP = indices_[j + 1];
     1091      value += duals[iRowM];
     1092      value -= duals[iRowP];
     1093      reducedCost[bestSequence] = value;
     1094      savedBestSequence_ = bestSequence;
     1095      savedBestDj_ = reducedCost[savedBestSequence_];
     1096    }
     1097  }
     1098  currentWanted_ = numberWanted;
    11161099}
    11171100// Allow any parts of a created CoinMatrix to be deleted
    1118 void
    1119 ClpNetworkMatrix::releasePackedMatrix() const
    1120 {
    1121      delete matrix_;
    1122      delete [] lengths_;
    1123      matrix_ = NULL;
    1124      lengths_ = NULL;
     1101void ClpNetworkMatrix::releasePackedMatrix() const
     1102{
     1103  delete matrix_;
     1104  delete[] lengths_;
     1105  matrix_ = NULL;
     1106  lengths_ = NULL;
    11251107}
    11261108// Append Columns
    1127 void
    1128 ClpNetworkMatrix::appendCols(int number, const CoinPackedVectorBase * const * columns)
    1129 {
    1130      int iColumn;
    1131      int numberBad = 0;
    1132      for (iColumn = 0; iColumn < number; iColumn++) {
    1133           int n = columns[iColumn]->getNumElements();
    1134           const double * element = columns[iColumn]->getElements();
    1135           if (n != 2)
    1136                numberBad++;
    1137           if (fabs(element[0]) != 1.0 || fabs(element[1]) != 1.0)
    1138                numberBad++;
    1139           else if (element[0]*element[1] != -1.0)
    1140                numberBad++;
    1141      }
    1142      if (numberBad)
    1143           throw CoinError("Not network", "appendCols", "ClpNetworkMatrix");
    1144      // Get rid of temporary arrays
    1145      delete [] lengths_;
    1146      lengths_ = NULL;
    1147      delete matrix_;
    1148      matrix_ = NULL;
    1149      CoinBigIndex size = 2 * number;
    1150      int * temp2 = new int [numberColumns_*2+size];
    1151      CoinMemcpyN(indices_, numberColumns_ * 2, temp2);
    1152      delete [] indices_;
    1153      indices_ = temp2;
    1154      // now add
    1155      size = 2 * numberColumns_;
    1156      for (iColumn = 0; iColumn < number; iColumn++) {
    1157           const int * row = columns[iColumn]->getIndices();
    1158           const double * element = columns[iColumn]->getElements();
    1159           if (element[0] == -1.0) {
    1160                indices_[size++] = row[0];
    1161                indices_[size++] = row[1];
    1162           } else {
    1163                indices_[size++] = row[1];
    1164                indices_[size++] = row[0];
    1165           }
    1166      }
    1167 
    1168      numberColumns_ += number;
     1109void ClpNetworkMatrix::appendCols(int number, const CoinPackedVectorBase *const *columns)
     1110{
     1111  int iColumn;
     1112  int numberBad = 0;
     1113  for (iColumn = 0; iColumn < number; iColumn++) {
     1114    int n = columns[iColumn]->getNumElements();
     1115    const double *element = columns[iColumn]->getElements();
     1116    if (n != 2)
     1117      numberBad++;
     1118    if (fabs(element[0]) != 1.0 || fabs(element[1]) != 1.0)
     1119      numberBad++;
     1120    else if (element[0] * element[1] != -1.0)
     1121      numberBad++;
     1122  }
     1123  if (numberBad)
     1124    throw CoinError("Not network", "appendCols", "ClpNetworkMatrix");
     1125  // Get rid of temporary arrays
     1126  delete[] lengths_;
     1127  lengths_ = NULL;
     1128  delete matrix_;
     1129  matrix_ = NULL;
     1130  CoinBigIndex size = 2 * number;
     1131  int *temp2 = new int[numberColumns_ * 2 + size];
     1132  CoinMemcpyN(indices_, numberColumns_ * 2, temp2);
     1133  delete[] indices_;
     1134  indices_ = temp2;
     1135  // now add
     1136  size = 2 * numberColumns_;
     1137  for (iColumn = 0; iColumn < number; iColumn++) {
     1138    const int *row = columns[iColumn]->getIndices();
     1139    const double *element = columns[iColumn]->getElements();
     1140    if (element[0] == -1.0) {
     1141      indices_[size++] = row[0];
     1142      indices_[size++] = row[1];
     1143    } else {
     1144      indices_[size++] = row[1];
     1145      indices_[size++] = row[0];
     1146    }
     1147  }
     1148
     1149  numberColumns_ += number;
    11691150}
    11701151// Append Rows
    1171 void
    1172 ClpNetworkMatrix::appendRows(int number, const CoinPackedVectorBase * const * rows)
    1173 {
    1174      // must be zero arrays
    1175      int numberBad = 0;
    1176      int iRow;
    1177      for (iRow = 0; iRow < number; iRow++) {
    1178           numberBad += rows[iRow]->getNumElements();
    1179      }
    1180      if (numberBad)
    1181           throw CoinError("Not NULL rows", "appendRows", "ClpNetworkMatrix");
    1182      numberRows_ += number;
     1152void ClpNetworkMatrix::appendRows(int number, const CoinPackedVectorBase *const *rows)
     1153{
     1154  // must be zero arrays
     1155  int numberBad = 0;
     1156  int iRow;
     1157  for (iRow = 0; iRow < number; iRow++) {
     1158    numberBad += rows[iRow]->getNumElements();
     1159  }
     1160  if (numberBad)
     1161    throw CoinError("Not NULL rows", "appendRows", "ClpNetworkMatrix");
     1162  numberRows_ += number;
    11831163}
    11841164#ifndef SLIM_CLP
     
    11871167   number of columns-1/rows-1 (if numberOther>0) or duplicates
    11881168   If 0 then rows, 1 if columns */
    1189 int
    1190 ClpNetworkMatrix::appendMatrix(int number, int type,
    1191                                const CoinBigIndex * starts, const int * index,
    1192                                const double * element, int /*numberOther*/)
    1193 {
    1194      int numberErrors = 0;
    1195      // make into CoinPackedVector
    1196      CoinPackedVectorBase ** vectors =
    1197           new CoinPackedVectorBase * [number];
    1198      int iVector;
    1199      for (iVector = 0; iVector < number; iVector++) {
    1200           CoinBigIndex iStart = starts[iVector];
    1201           vectors[iVector] =
    1202             new CoinPackedVector(static_cast<int>(starts[iVector+1] - iStart),
    1203                                     index + iStart, element + iStart);
    1204      }
    1205      if (type == 0) {
    1206           // rows
    1207           appendRows(number, vectors);
    1208      } else {
    1209           // columns
    1210           appendCols(number, vectors);
    1211      }
    1212      for (iVector = 0; iVector < number; iVector++)
    1213           delete vectors[iVector];
    1214      delete [] vectors;
    1215      return numberErrors;
     1169int ClpNetworkMatrix::appendMatrix(int number, int type,
     1170  const CoinBigIndex *starts, const int *index,
     1171  const double *element, int /*numberOther*/)
     1172{
     1173  int numberErrors = 0;
     1174  // make into CoinPackedVector
     1175  CoinPackedVectorBase **vectors = new CoinPackedVectorBase *[number];
     1176  int iVector;
     1177  for (iVector = 0; iVector < number; iVector++) {
     1178    CoinBigIndex iStart = starts[iVector];
     1179    vectors[iVector] = new CoinPackedVector(static_cast< int >(starts[iVector + 1] - iStart),
     1180      index + iStart, element + iStart);
     1181  }
     1182  if (type == 0) {
     1183    // rows
     1184    appendRows(number, vectors);
     1185  } else {
     1186    // columns
     1187    appendCols(number, vectors);
     1188  }
     1189  for (iVector = 0; iVector < number; iVector++)
     1190    delete vectors[iVector];
     1191  delete[] vectors;
     1192  return numberErrors;
    12161193}
    12171194#endif
     
    12191196   and order is as given */
    12201197ClpMatrixBase *
    1221 ClpNetworkMatrix::subsetClone (int numberRows, const int * whichRows,
    1222                                int numberColumns,
    1223                                const int * whichColumns) const
    1224 {
    1225      return new ClpNetworkMatrix(*this, numberRows, whichRows,
    1226                                  numberColumns, whichColumns);
     1198ClpNetworkMatrix::subsetClone(int numberRows, const int *whichRows,
     1199  int numberColumns,
     1200  const int *whichColumns) const
     1201{
     1202  return new ClpNetworkMatrix(*this, numberRows, whichRows,
     1203    numberColumns, whichColumns);
    12271204}
    12281205/* Subset constructor (without gaps).  Duplicates are allowed
    12291206   and order is as given */
    1230 ClpNetworkMatrix::ClpNetworkMatrix (
    1231      const ClpNetworkMatrix & rhs,
    1232      int numberRows, const int * whichRow,
    1233      int numberColumns, const int * whichColumn)
    1234      : ClpMatrixBase(rhs)
    1235 {
    1236      setType(11);
    1237      matrix_ = NULL;
    1238      lengths_ = NULL;
    1239      indices_ = new int[2*numberColumns];;
    1240      numberRows_ = numberRows;
    1241      numberColumns_ = numberColumns;
    1242      trueNetwork_ = true;
    1243      int iColumn;
    1244      int numberBad = 0;
    1245      int * which = new int [rhs.numberRows_];
    1246      int iRow;
    1247      for (iRow = 0; iRow < rhs.numberRows_; iRow++)
    1248           which[iRow] = -1;
    1249      int n = 0;
    1250      for (iRow = 0; iRow < numberRows; iRow++) {
    1251           int jRow = whichRow[iRow];
    1252           assert (jRow >= 0 && jRow < rhs.numberRows_);
    1253           which[jRow] = n++;
    1254      }
    1255      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    1256           CoinBigIndex start;
    1257           CoinBigIndex i;
    1258           start = 2 * iColumn;
    1259           CoinBigIndex offset = 2 * whichColumn[iColumn] - start;
    1260           for (i = start; i < start + 2; i++) {
    1261                int iRow = rhs.indices_[i+offset];
    1262                iRow = which[iRow];
    1263                if (iRow < 0)
    1264                     numberBad++;
    1265                else
    1266                     indices_[i] = iRow;
    1267           }
    1268      }
    1269      if (numberBad)
    1270           throw CoinError("Invalid rows", "subsetConstructor", "ClpNetworkMatrix");
    1271 }
     1207ClpNetworkMatrix::ClpNetworkMatrix(
     1208  const ClpNetworkMatrix &rhs,
     1209  int numberRows, const int *whichRow,
     1210  int numberColumns, const int *whichColumn)
     1211  : ClpMatrixBase(rhs)
     1212{
     1213  setType(11);
     1214  matrix_ = NULL;
     1215  lengths_ = NULL;
     1216  indices_ = new int[2 * numberColumns];
     1217  ;
     1218  numberRows_ = numberRows;
     1219  numberColumns_ = numberColumns;
     1220  trueNetwork_ = true;
     1221  int iColumn;
     1222  int numberBad = 0;
     1223  int *which = new int[rhs.numberRows_];
     1224  int iRow;
     1225  for (iRow = 0; iRow < rhs.numberRows_; iRow++)
     1226    which[iRow] = -1;
     1227  int n = 0;
     1228  for (iRow = 0; iRow < numberRows; iRow++) {
     1229    int jRow = whichRow[iRow];
     1230    assert(jRow >= 0 && jRow < rhs.numberRows_);
     1231    which[jRow] = n++;
     1232  }
     1233  for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     1234    CoinBigIndex start;
     1235    CoinBigIndex i;
     1236    start = 2 * iColumn;
     1237    CoinBigIndex offset = 2 * whichColumn[iColumn] - start;
     1238    for (i = start; i < start + 2; i++) {
     1239      int iRow = rhs.indices_[i + offset];
     1240      iRow = which[iRow];
     1241      if (iRow < 0)
     1242        numberBad++;
     1243      else
     1244        indices_[i] = iRow;
     1245    }
     1246  }
     1247  if (numberBad)
     1248    throw CoinError("Invalid rows", "subsetConstructor", "ClpNetworkMatrix");
     1249}
     1250
     1251/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
     1252*/
Note: See TracChangeset for help on using the changeset viewer.