Ignore:
Timestamp:
Feb 26, 2010 12:27:59 PM (10 years ago)
Author:
mjs
Message:

Formatted .cpp, .hpp, .c, .h files with "astyle -A4 -p". This matches the formatting used in the grand CBC reorganization.

File:
1 edited

Legend:

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

    r1402 r1525  
    2525
    2626//-------------------------------------------------------------------
    27 // Default Constructor 
     27// Default Constructor
    2828//-------------------------------------------------------------------
    29 ClpNetworkMatrix::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;
     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;
    3939}
    4040
    4141/* Constructor from two arrays */
    4242ClpNetworkMatrix::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   numberRows_=-1;
    51   numberColumns_=numberColumns;
    52   trueNetwork_=true;
    53   int iColumn;
    54   CoinBigIndex j=0;
    55   for (iColumn=0;iColumn<numberColumns_;iColumn++, j+=2) {
    56     int iRow = head[iColumn];
    57     numberRows_ = CoinMax(numberRows_,iRow);
    58     indices_[j]=iRow;
    59     iRow = tail[iColumn];
    60     numberRows_ = CoinMax(numberRows_,iRow);
    61     indices_[j+1]=iRow;
    62   }
    63   numberRows_++;
     43                                   const int * tail)
     44     : ClpMatrixBase()
     45{
     46     setType(11);
     47     matrix_ = NULL;
     48     lengths_ = NULL;
     49     indices_ = new int[2*numberColumns];;
     50     numberRows_ = -1;
     51     numberColumns_ = numberColumns;
     52     trueNetwork_ = true;
     53     int iColumn;
     54     CoinBigIndex j = 0;
     55     for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
     56          int iRow = head[iColumn];
     57          numberRows_ = CoinMax(numberRows_, iRow);
     58          indices_[j] = iRow;
     59          iRow = tail[iColumn];
     60          numberRows_ = CoinMax(numberRows_, iRow);
     61          indices_[j+1] = iRow;
     62     }
     63     numberRows_++;
    6464}
    6565//-------------------------------------------------------------------
    66 // Copy constructor 
     66// Copy constructor
    6767//-------------------------------------------------------------------
    68 ClpNetworkMatrix::ClpNetworkMatrix (const ClpNetworkMatrix & rhs) 
    69 : ClpMatrixBase(rhs)
    70 { 
    71   matrix_ = NULL;
    72   lengths_=NULL;
    73   indices_=NULL;
    74   numberRows_=rhs.numberRows_;
    75   numberColumns_=rhs.numberColumns_;
    76   trueNetwork_=rhs.trueNetwork_;
    77   if (numberColumns_) {
    78     indices_ = new int [ 2*numberColumns_];
    79     CoinMemcpyN(rhs.indices_,2*numberColumns_,indices_);
    80   }
    81   int numberRows = getNumRows();
    82   if (rhs.rhsOffset_&&numberRows) {
    83     rhsOffset_ = ClpCopyOfArray(rhs.rhsOffset_,numberRows);
    84   } else {
    85     rhsOffset_=NULL;
    86   }
    87 }
    88 
    89 ClpNetworkMatrix::ClpNetworkMatrix (const CoinPackedMatrix & rhs) 
    90   : ClpMatrixBase()
    91 { 
    92   setType(11);
    93   matrix_ = NULL;
    94   lengths_=NULL;
    95   indices_=NULL;
    96   int iColumn;
    97   assert (rhs.isColOrdered());
    98   // get matrix data pointers
    99   const int * row = rhs.getIndices();
    100   const CoinBigIndex * columnStart = rhs.getVectorStarts();
    101   const int * columnLength = rhs.getVectorLengths();
    102   const double * elementByColumn = rhs.getElements();
    103   numberColumns_ = rhs.getNumCols();
    104   int goodNetwork=1;
    105   numberRows_=-1;
    106   indices_ = new int[2*numberColumns_];
    107   CoinBigIndex j=0;
    108   for (iColumn=0;iColumn<numberColumns_;iColumn++, j+=2) {
    109     CoinBigIndex k=columnStart[iColumn];
    110     int iRow;
    111     switch (columnLength[iColumn]) {
    112     case 0:
    113       goodNetwork=-1; // not classic network
    114       indices_[j]=-1;
    115       indices_[j+1]=-1;
    116       break;
    117      
    118     case 1:
    119       goodNetwork=-1; // not classic network
    120       if (fabs(elementByColumn[k]-1.0)<1.0e-10) {
    121         indices_[j] = -1;
    122         iRow = row[k];
    123         numberRows_ = CoinMax(numberRows_,iRow);
    124         indices_[j+1]=iRow;
    125       } else if (fabs(elementByColumn[k]+1.0)<1.0e-10) {
    126         indices_[j+1] = -1;
    127         iRow = row[k];
    128         numberRows_ = CoinMax(numberRows_,iRow);
    129         indices_[j]=iRow;
    130       } else {
    131         goodNetwork = 0; // not a network
    132       }
    133       break;
    134      
    135     case 2:
    136       if (fabs(elementByColumn[k]-1.0)<1.0e-10) {
    137         if (fabs(elementByColumn[k+1]+1.0)<1.0e-10) {
    138           iRow = row[k];
    139           numberRows_ = CoinMax(numberRows_,iRow);
    140           indices_[j+1]=iRow;
    141           iRow = row[k+1];
    142           numberRows_ = CoinMax(numberRows_,iRow);
    143           indices_[j] = iRow;
    144         } else {
    145           goodNetwork = 0; // not a network
    146         }
    147       } else if (fabs(elementByColumn[k]+1.0)<1.0e-10) {
    148         if (fabs(elementByColumn[k+1]-1.0)<1.0e-10) {
    149           iRow = row[k];
    150           numberRows_ = CoinMax(numberRows_,iRow);
    151           indices_[j]=iRow;
    152           iRow = row[k+1];
    153           numberRows_ = CoinMax(numberRows_,iRow);
    154           indices_[j+1] = iRow;
    155         } else {
    156           goodNetwork = 0; // not a network
    157         }
    158       } else {
    159         goodNetwork = 0; // not a network
    160       }
    161       break;
    162 
    163     default:
    164       goodNetwork = 0; // not a network
    165       break;
    166     }
    167     if (!goodNetwork)
    168       break;
    169   }
    170   if (!goodNetwork) {
    171     delete [] indices_;
    172     // put in message
    173     printf("Not a network - can test if indices_ null\n");
    174     indices_=NULL;
    175     numberRows_=0;
    176     numberColumns_=0;
    177   } else {
    178     numberRows_ ++; //  correct
    179     trueNetwork_ = goodNetwork>0;
    180   }
     68ClpNetworkMatrix::ClpNetworkMatrix (const ClpNetworkMatrix & rhs)
     69     : ClpMatrixBase(rhs)
     70{
     71     matrix_ = NULL;
     72     lengths_ = NULL;
     73     indices_ = NULL;
     74     numberRows_ = rhs.numberRows_;
     75     numberColumns_ = rhs.numberColumns_;
     76     trueNetwork_ = rhs.trueNetwork_;
     77     if (numberColumns_) {
     78          indices_ = new int [ 2*numberColumns_];
     79          CoinMemcpyN(rhs.indices_, 2 * numberColumns_, indices_);
     80     }
     81     int numberRows = getNumRows();
     82     if (rhs.rhsOffset_ && numberRows) {
     83          rhsOffset_ = ClpCopyOfArray(rhs.rhsOffset_, numberRows);
     84     } else {
     85          rhsOffset_ = NULL;
     86     }
     87}
     88
     89ClpNetworkMatrix::ClpNetworkMatrix (const CoinPackedMatrix & rhs)
     90     : ClpMatrixBase()
     91{
     92     setType(11);
     93     matrix_ = NULL;
     94     lengths_ = NULL;
     95     indices_ = NULL;
     96     int iColumn;
     97     assert (rhs.isColOrdered());
     98     // get matrix data pointers
     99     const int * row = rhs.getIndices();
     100     const CoinBigIndex * columnStart = rhs.getVectorStarts();
     101     const int * columnLength = rhs.getVectorLengths();
     102     const double * elementByColumn = rhs.getElements();
     103     numberColumns_ = rhs.getNumCols();
     104     int goodNetwork = 1;
     105     numberRows_ = -1;
     106     indices_ = new int[2*numberColumns_];
     107     CoinBigIndex j = 0;
     108     for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
     109          CoinBigIndex k = columnStart[iColumn];
     110          int iRow;
     111          switch (columnLength[iColumn]) {
     112          case 0:
     113               goodNetwork = -1; // not classic network
     114               indices_[j] = -1;
     115               indices_[j+1] = -1;
     116               break;
     117
     118          case 1:
     119               goodNetwork = -1; // not classic network
     120               if (fabs(elementByColumn[k] - 1.0) < 1.0e-10) {
     121                    indices_[j] = -1;
     122                    iRow = row[k];
     123                    numberRows_ = CoinMax(numberRows_, iRow);
     124                    indices_[j+1] = iRow;
     125               } else if (fabs(elementByColumn[k] + 1.0) < 1.0e-10) {
     126                    indices_[j+1] = -1;
     127                    iRow = row[k];
     128                    numberRows_ = CoinMax(numberRows_, iRow);
     129                    indices_[j] = iRow;
     130               } else {
     131                    goodNetwork = 0; // not a network
     132               }
     133               break;
     134
     135          case 2:
     136               if (fabs(elementByColumn[k] - 1.0) < 1.0e-10) {
     137                    if (fabs(elementByColumn[k+1] + 1.0) < 1.0e-10) {
     138                         iRow = row[k];
     139                         numberRows_ = CoinMax(numberRows_, iRow);
     140                         indices_[j+1] = iRow;
     141                         iRow = row[k+1];
     142                         numberRows_ = CoinMax(numberRows_, iRow);
     143                         indices_[j] = iRow;
     144                    } else {
     145                         goodNetwork = 0; // not a network
     146                    }
     147               } else if (fabs(elementByColumn[k] + 1.0) < 1.0e-10) {
     148                    if (fabs(elementByColumn[k+1] - 1.0) < 1.0e-10) {
     149                         iRow = row[k];
     150                         numberRows_ = CoinMax(numberRows_, iRow);
     151                         indices_[j] = iRow;
     152                         iRow = row[k+1];
     153                         numberRows_ = CoinMax(numberRows_, iRow);
     154                         indices_[j+1] = iRow;
     155                    } else {
     156                         goodNetwork = 0; // not a network
     157                    }
     158               } else {
     159                    goodNetwork = 0; // not a network
     160               }
     161               break;
     162
     163          default:
     164               goodNetwork = 0; // not a network
     165               break;
     166          }
     167          if (!goodNetwork)
     168               break;
     169     }
     170     if (!goodNetwork) {
     171          delete [] indices_;
     172          // put in message
     173          printf("Not a network - can test if indices_ null\n");
     174          indices_ = NULL;
     175          numberRows_ = 0;
     176          numberColumns_ = 0;
     177     } else {
     178          numberRows_ ++; //  correct
     179          trueNetwork_ = goodNetwork > 0;
     180     }
    181181}
    182182
    183183//-------------------------------------------------------------------
    184 // Destructor 
     184// Destructor
    185185//-------------------------------------------------------------------
    186186ClpNetworkMatrix::~ClpNetworkMatrix ()
    187187{
    188   delete matrix_;
    189   delete [] lengths_;
    190   delete [] indices_;
     188     delete matrix_;
     189     delete [] lengths_;
     190     delete [] indices_;
    191191}
    192192
    193193//----------------------------------------------------------------
    194 // Assignment operator 
     194// Assignment operator
    195195//-------------------------------------------------------------------
    196196ClpNetworkMatrix &
    197197ClpNetworkMatrix::operator=(const ClpNetworkMatrix& rhs)
    198198{
    199   if (this != &rhs) {
    200     ClpMatrixBase::operator=(rhs);
    201     delete matrix_;
    202     delete [] lengths_;
    203     delete [] indices_;
    204     matrix_ = NULL;
    205     lengths_=NULL;
    206     indices_=NULL;
    207     numberRows_=rhs.numberRows_;
    208     numberColumns_=rhs.numberColumns_;
    209     trueNetwork_=rhs.trueNetwork_;
    210     if (numberColumns_) {
    211       indices_ = new int [ 2*numberColumns_];
    212       CoinMemcpyN(rhs.indices_,2*numberColumns_,indices_);
    213     }
    214   }
    215   return *this;
     199     if (this != &rhs) {
     200          ClpMatrixBase::operator=(rhs);
     201          delete matrix_;
     202          delete [] lengths_;
     203          delete [] indices_;
     204          matrix_ = NULL;
     205          lengths_ = NULL;
     206          indices_ = NULL;
     207          numberRows_ = rhs.numberRows_;
     208          numberColumns_ = rhs.numberColumns_;
     209          trueNetwork_ = rhs.trueNetwork_;
     210          if (numberColumns_) {
     211               indices_ = new int [ 2*numberColumns_];
     212               CoinMemcpyN(rhs.indices_, 2 * numberColumns_, indices_);
     213          }
     214     }
     215     return *this;
    216216}
    217217//-------------------------------------------------------------------
     
    220220ClpMatrixBase * ClpNetworkMatrix::clone() const
    221221{
    222   return new ClpNetworkMatrix(*this);
     222     return new ClpNetworkMatrix(*this);
    223223}
    224224
    225225/* Returns a new matrix in reverse order without gaps */
    226 ClpMatrixBase * 
     226ClpMatrixBase *
    227227ClpNetworkMatrix::reverseOrderedCopy() const
    228228{
    229   // count number in each row
    230   CoinBigIndex * tempP = new CoinBigIndex [numberRows_];
    231   CoinBigIndex * tempN = new CoinBigIndex [numberRows_];
    232   memset(tempP,0,numberRows_*sizeof(CoinBigIndex));
    233   memset(tempN,0,numberRows_*sizeof(CoinBigIndex));
    234   CoinBigIndex j=0;
    235   int i;
    236   for (i=0;i<numberColumns_;i++,j+=2) {
    237     int iRow = indices_[j];
    238     tempN[iRow]++;
    239     iRow = indices_[j+1];
    240     tempP[iRow]++;
    241   }
    242   int * newIndices = new int [2*numberColumns_];
    243   CoinBigIndex * newP = new CoinBigIndex [numberRows_+1];
    244   CoinBigIndex * newN = new CoinBigIndex[numberRows_];
    245   int iRow;
    246   j=0;
    247   // do starts
    248   for (iRow=0;iRow<numberRows_;iRow++) {
    249     newP[iRow]=j;
    250     j += tempP[iRow];
    251     tempP[iRow]=newP[iRow];
    252     newN[iRow] = j;
    253     j += tempN[iRow];
    254     tempN[iRow]=newN[iRow];
    255   }
    256   newP[numberRows_]=j;
    257   j=0;
    258   for (i=0;i<numberColumns_;i++,j+=2) {
    259     int iRow = indices_[j];
    260     CoinBigIndex put = tempN[iRow];
    261     newIndices[put++] = i;
    262     tempN[iRow] = put;
    263     iRow = indices_[j+1];
    264     put = tempP[iRow];
    265     newIndices[put++] = i;
    266     tempP[iRow] = put;
    267   }
    268   delete [] tempP;
    269   delete [] tempN;
    270   ClpPlusMinusOneMatrix * newCopy = new ClpPlusMinusOneMatrix();
    271   newCopy->passInCopy(numberRows_, numberColumns_,
    272                       false,  newIndices, newP, newN);
    273   return newCopy;
     229     // count number in each row
     230     CoinBigIndex * tempP = new CoinBigIndex [numberRows_];
     231     CoinBigIndex * tempN = new CoinBigIndex [numberRows_];
     232     memset(tempP, 0, numberRows_ * sizeof(CoinBigIndex));
     233     memset(tempN, 0, numberRows_ * sizeof(CoinBigIndex));
     234     CoinBigIndex j = 0;
     235     int i;
     236     for (i = 0; i < numberColumns_; i++, j += 2) {
     237          int iRow = indices_[j];
     238          tempN[iRow]++;
     239          iRow = indices_[j+1];
     240          tempP[iRow]++;
     241     }
     242     int * newIndices = new int [2*numberColumns_];
     243     CoinBigIndex * newP = new CoinBigIndex [numberRows_+1];
     244     CoinBigIndex * newN = new CoinBigIndex[numberRows_];
     245     int iRow;
     246     j = 0;
     247     // do starts
     248     for (iRow = 0; iRow < numberRows_; iRow++) {
     249          newP[iRow] = j;
     250          j += tempP[iRow];
     251          tempP[iRow] = newP[iRow];
     252          newN[iRow] = j;
     253          j += tempN[iRow];
     254          tempN[iRow] = newN[iRow];
     255     }
     256     newP[numberRows_] = j;
     257     j = 0;
     258     for (i = 0; i < numberColumns_; i++, j += 2) {
     259          int iRow = indices_[j];
     260          CoinBigIndex put = tempN[iRow];
     261          newIndices[put++] = i;
     262          tempN[iRow] = put;
     263          iRow = indices_[j+1];
     264          put = tempP[iRow];
     265          newIndices[put++] = i;
     266          tempP[iRow] = put;
     267     }
     268     delete [] tempP;
     269     delete [] tempN;
     270     ClpPlusMinusOneMatrix * newCopy = new ClpPlusMinusOneMatrix();
     271     newCopy->passInCopy(numberRows_, numberColumns_,
     272                         false,  newIndices, newP, newN);
     273     return newCopy;
    274274}
    275275//unscaled versions
    276 void 
     276void
    277277ClpNetworkMatrix::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 }
    307 void 
     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
    308308ClpNetworkMatrix::transposeTimes(double scalar,
    309                                 const double * x, double * y) const
    310 {
    311   int iColumn;
    312   CoinBigIndex j=0;
    313   if (trueNetwork_) {
    314     for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
    315       double value = y[iColumn];
    316       int iRowM = indices_[j];
    317       int iRowP = indices_[j+1];
    318       value -= scalar*x[iRowM];
    319       value += scalar*x[iRowP];
    320       y[iColumn] = value;
    321     }
    322   } else {
    323     // skip negative rows
    324     for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
    325       double value = y[iColumn];
    326       int iRowM = indices_[j];
    327       int iRowP = indices_[j+1];
    328       if (iRowM>=0)
    329         value -= scalar*x[iRowM];
    330       if (iRowP>=0)
    331         value += scalar*x[iRowP];
    332       y[iColumn] = value;
    333     }
    334   }
    335 }
    336 void 
     309                                 const double * x, double * y) const
     310{
     311     int iColumn;
     312     CoinBigIndex j = 0;
     313     if (trueNetwork_) {
     314          for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
     315               double value = y[iColumn];
     316               int iRowM = indices_[j];
     317               int iRowP = indices_[j+1];
     318               value -= scalar * x[iRowM];
     319               value += scalar * x[iRowP];
     320               y[iColumn] = value;
     321          }
     322     } else {
     323          // skip negative rows
     324          for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
     325               double value = y[iColumn];
     326               int iRowM = indices_[j];
     327               int iRowP = indices_[j+1];
     328               if (iRowM >= 0)
     329                    value -= scalar * x[iRowM];
     330               if (iRowP >= 0)
     331                    value += scalar * x[iRowP];
     332               y[iColumn] = value;
     333          }
     334     }
     335}
     336void
    337337ClpNetworkMatrix::times(double scalar,
    338                        const double * x, double * y,
    339                         const double * /*rowScale*/,
    340                         const double * /*columnScale*/) const
    341 {
    342   // we know it is not scaled
    343   times(scalar, x, y);
    344 }
    345 void 
     338                        const double * x, double * y,
     339                        const double * /*rowScale*/,
     340                        const double * /*columnScale*/) const
     341{
     342     // we know it is not scaled
     343     times(scalar, x, y);
     344}
     345void
    346346ClpNetworkMatrix::transposeTimes( double scalar,
    347                                 const double * x, double * y,
    348                                   const double * /*rowScale*/,
    349                                   const double * /*columnScale*/,
    350                                   double * /*spare*/) const
    351 {
    352   // we know it is not scaled
    353   transposeTimes(scalar, x, y);
    354 }
    355 /* Return <code>x * A + y</code> in <code>z</code>. 
     347                                  const double * x, double * y,
     348                                  const double * /*rowScale*/,
     349                                  const double * /*columnScale*/,
     350                                  double * /*spare*/) const
     351{
     352     // we know it is not scaled
     353     transposeTimes(scalar, x, y);
     354}
     355/* Return <code>x * A + y</code> in <code>z</code>.
    356356        Squashes small elements and knows about ClpSimplex */
    357 void 
     357void
    358358ClpNetworkMatrix::transposeTimes(const ClpSimplex * model, double scalar,
    359                               const CoinIndexedVector * rowArray,
    360                               CoinIndexedVector * y,
    361                               CoinIndexedVector * columnArray) const
    362 {
    363   // we know it is not scaled
    364   columnArray->clear();
    365   double * pi = rowArray->denseVector();
    366   int numberNonZero=0;
    367   int * index = columnArray->getIndices();
    368   double * array = columnArray->denseVector();
    369   int numberInRowArray = rowArray->getNumElements();
    370   // maybe I need one in OsiSimplex
    371   double zeroTolerance = model->zeroTolerance();
    372   int numberRows = model->numberRows();
     359                                 const CoinIndexedVector * rowArray,
     360                                 CoinIndexedVector * y,
     361                                 CoinIndexedVector * columnArray) const
     362{
     363     // we know it is not scaled
     364     columnArray->clear();
     365     double * pi = rowArray->denseVector();
     366     int numberNonZero = 0;
     367     int * index = columnArray->getIndices();
     368     double * array = columnArray->denseVector();
     369     int numberInRowArray = rowArray->getNumElements();
     370     // maybe I need one in OsiSimplex
     371     double zeroTolerance = model->zeroTolerance();
     372     int numberRows = model->numberRows();
    373373#ifndef NO_RTTI
    374   ClpPlusMinusOneMatrix* rowCopy =
    375     dynamic_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
     374     ClpPlusMinusOneMatrix* rowCopy =
     375          dynamic_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
    376376#else
    377   ClpPlusMinusOneMatrix* rowCopy =
    378     static_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
     377     ClpPlusMinusOneMatrix* rowCopy =
     378          static_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
    379379#endif
    380   bool packed = rowArray->packedMode();
    381   double factor = 0.3;
    382   // We may not want to do by row if there may be cache problems
    383   int numberColumns = model->numberColumns();
    384   // It would be nice to find L2 cache size - for moment 512K
    385   // Be slightly optimistic
    386   if (numberColumns*sizeof(double)>1000000) {
    387     if (numberRows*10<numberColumns)
    388       factor=0.1;
    389     else if (numberRows*4<numberColumns)
    390       factor=0.15;
    391     else if (numberRows*2<numberColumns)
    392       factor=0.2;
    393     //if (model->numberIterations()%50==0)
    394     //printf("%d nonzero\n",numberInRowArray);
    395   }
    396   if (numberInRowArray>factor*numberRows||!rowCopy) {
    397     // do by column
    398     int iColumn;
    399     assert (!y->getNumElements());
    400     CoinBigIndex j=0;
    401     if (packed) {
    402       // need to expand pi into y
    403       assert(y->capacity()>=numberRows);
    404       double * piOld = pi;
    405       pi = y->denseVector();
    406       const int * whichRow = rowArray->getIndices();
    407       int i;
    408       // modify pi so can collapse to one loop
    409       for (i=0;i<numberInRowArray;i++) {
    410         int iRow = whichRow[i];
    411         pi[iRow]=scalar*piOld[i];
    412       }
    413       if (trueNetwork_) {
    414         for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
    415           double value = 0.0;
    416           int iRowM = indices_[j];
    417           int iRowP = indices_[j+1];
    418           value -= pi[iRowM];
    419           value += pi[iRowP];
    420           if (fabs(value)>zeroTolerance) {
    421             array[numberNonZero]=value;
    422             index[numberNonZero++]=iColumn;
    423           }
    424         }
    425       } else {
    426         // skip negative rows
    427         for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
    428           double value = 0.0;
    429           int iRowM = indices_[j];
    430           int iRowP = indices_[j+1];
    431           if (iRowM>=0)
    432             value -= pi[iRowM];
    433           if (iRowP>=0)
    434             value += pi[iRowP];
    435           if (fabs(value)>zeroTolerance) {
    436             array[numberNonZero]=value;
    437             index[numberNonZero++]=iColumn;
    438           }
    439         }
    440       }
    441       for (i=0;i<numberInRowArray;i++) {
    442         int iRow = whichRow[i];
    443         pi[iRow]=0.0;
    444       }
    445     } else {
    446       if (trueNetwork_) {
    447         for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
    448           double value = 0.0;
    449           int iRowM = indices_[j];
    450           int iRowP = indices_[j+1];
    451           value -= scalar*pi[iRowM];
    452           value += scalar*pi[iRowP];
    453           if (fabs(value)>zeroTolerance) {
    454             index[numberNonZero++]=iColumn;
    455             array[iColumn]=value;
    456           }
    457         }
    458       } else {
    459         // skip negative rows
    460         for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
    461           double value = 0.0;
    462           int iRowM = indices_[j];
    463           int iRowP = indices_[j+1];
    464           if (iRowM>=0)
    465             value -= scalar*pi[iRowM];
    466           if (iRowP>=0)
    467             value += scalar*pi[iRowP];
    468           if (fabs(value)>zeroTolerance) {
    469             index[numberNonZero++]=iColumn;
    470             array[iColumn]=value;
    471           }
    472         }
    473       }
    474     }
    475     columnArray->setNumElements(numberNonZero);
    476   } else {
    477     // do by row
    478     rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray);
    479   }
     380     bool packed = rowArray->packedMode();
     381     double factor = 0.3;
     382     // We may not want to do by row if there may be cache problems
     383     int numberColumns = model->numberColumns();
     384     // It would be nice to find L2 cache size - for moment 512K
     385     // Be slightly optimistic
     386     if (numberColumns * sizeof(double) > 1000000) {
     387          if (numberRows * 10 < numberColumns)
     388               factor = 0.1;
     389          else if (numberRows * 4 < numberColumns)
     390               factor = 0.15;
     391          else if (numberRows * 2 < numberColumns)
     392               factor = 0.2;
     393          //if (model->numberIterations()%50==0)
     394          //printf("%d nonzero\n",numberInRowArray);
     395     }
     396     if (numberInRowArray > factor * numberRows || !rowCopy) {
     397          // do by column
     398          int iColumn;
     399          assert (!y->getNumElements());
     400          CoinBigIndex j = 0;
     401          if (packed) {
     402               // need to expand pi into y
     403               assert(y->capacity() >= numberRows);
     404               double * piOld = pi;
     405               pi = y->denseVector();
     406               const int * whichRow = rowArray->getIndices();
     407               int i;
     408               // modify pi so can collapse to one loop
     409               for (i = 0; i < numberInRowArray; i++) {
     410                    int iRow = whichRow[i];
     411                    pi[iRow] = scalar * piOld[i];
     412               }
     413               if (trueNetwork_) {
     414                    for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
     415                         double value = 0.0;
     416                         int iRowM = indices_[j];
     417                         int iRowP = indices_[j+1];
     418                         value -= pi[iRowM];
     419                         value += pi[iRowP];
     420                         if (fabs(value) > zeroTolerance) {
     421                              array[numberNonZero] = value;
     422                              index[numberNonZero++] = iColumn;
     423                         }
     424                    }
     425               } else {
     426                    // skip negative rows
     427                    for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
     428                         double value = 0.0;
     429                         int iRowM = indices_[j];
     430                         int iRowP = indices_[j+1];
     431                         if (iRowM >= 0)
     432                              value -= pi[iRowM];
     433                         if (iRowP >= 0)
     434                              value += pi[iRowP];
     435                         if (fabs(value) > zeroTolerance) {
     436                              array[numberNonZero] = value;
     437                              index[numberNonZero++] = iColumn;
     438                         }
     439                    }
     440               }
     441               for (i = 0; i < numberInRowArray; i++) {
     442                    int iRow = whichRow[i];
     443                    pi[iRow] = 0.0;
     444               }
     445          } else {
     446               if (trueNetwork_) {
     447                    for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
     448                         double value = 0.0;
     449                         int iRowM = indices_[j];
     450                         int iRowP = indices_[j+1];
     451                         value -= scalar * pi[iRowM];
     452                         value += scalar * pi[iRowP];
     453                         if (fabs(value) > zeroTolerance) {
     454                              index[numberNonZero++] = iColumn;
     455                              array[iColumn] = value;
     456                         }
     457                    }
     458               } else {
     459                    // skip negative rows
     460                    for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
     461                         double value = 0.0;
     462                         int iRowM = indices_[j];
     463                         int iRowP = indices_[j+1];
     464                         if (iRowM >= 0)
     465                              value -= scalar * pi[iRowM];
     466                         if (iRowP >= 0)
     467                              value += scalar * pi[iRowP];
     468                         if (fabs(value) > zeroTolerance) {
     469                              index[numberNonZero++] = iColumn;
     470                              array[iColumn] = value;
     471                         }
     472                    }
     473               }
     474          }
     475          columnArray->setNumElements(numberNonZero);
     476     } else {
     477          // do by row
     478          rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray);
     479     }
    480480}
    481481/* Return <code>x *A in <code>z</code> but
    482482   just for indices in y. */
    483 void 
     483void
    484484ClpNetworkMatrix::subsetTransposeTimes(const ClpSimplex * /*model*/,
    485                               const CoinIndexedVector * rowArray,
    486                               const CoinIndexedVector * y,
    487                               CoinIndexedVector * columnArray) const
    488 {
    489   columnArray->clear();
    490   double * pi = rowArray->denseVector();
    491   double * array = columnArray->denseVector();
    492   int jColumn;
    493   int numberToDo = y->getNumElements();
    494   const int * which = y->getIndices();
    495   assert (!rowArray->packedMode());
    496   columnArray->setPacked();
    497   if (trueNetwork_) {
    498     for (jColumn=0;jColumn<numberToDo;jColumn++) {
    499       int iColumn = which[jColumn];
    500       double value = 0.0;
    501       CoinBigIndex j=iColumn<<1;
    502       int iRowM = indices_[j];
    503       int iRowP = indices_[j+1];
    504       value -= pi[iRowM];
    505       value += pi[iRowP];
    506       array[jColumn]=value;
    507     }
    508   } else {
    509     // skip negative rows
    510     for (jColumn=0;jColumn<numberToDo;jColumn++) {
    511       int iColumn = which[jColumn];
    512       double value = 0.0;
    513       CoinBigIndex j=iColumn<<1;
    514       int iRowM = indices_[j];
    515       int iRowP = indices_[j+1];
    516       if (iRowM>=0)
    517         value -= pi[iRowM];
    518       if (iRowP>=0)
    519         value += pi[iRowP];
    520       array[jColumn]=value;
    521     }
    522   }
     485                                       const CoinIndexedVector * rowArray,
     486                                       const CoinIndexedVector * y,
     487                                       CoinIndexedVector * columnArray) const
     488{
     489     columnArray->clear();
     490     double * pi = rowArray->denseVector();
     491     double * array = columnArray->denseVector();
     492     int jColumn;
     493     int numberToDo = y->getNumElements();
     494     const int * which = y->getIndices();
     495     assert (!rowArray->packedMode());
     496     columnArray->setPacked();
     497     if (trueNetwork_) {
     498          for (jColumn = 0; jColumn < numberToDo; jColumn++) {
     499               int iColumn = which[jColumn];
     500               double value = 0.0;
     501               CoinBigIndex j = iColumn << 1;
     502               int iRowM = indices_[j];
     503               int iRowP = indices_[j+1];
     504               value -= pi[iRowM];
     505               value += pi[iRowP];
     506               array[jColumn] = value;
     507          }
     508     } else {
     509          // skip negative rows
     510          for (jColumn = 0; jColumn < numberToDo; jColumn++) {
     511               int iColumn = which[jColumn];
     512               double value = 0.0;
     513               CoinBigIndex j = iColumn << 1;
     514               int iRowM = indices_[j];
     515               int iRowP = indices_[j+1];
     516               if (iRowM >= 0)
     517                    value -= pi[iRowM];
     518               if (iRowP >= 0)
     519                    value += pi[iRowP];
     520               array[jColumn] = value;
     521          }
     522     }
    523523}
    524524/// returns number of elements in column part of basis,
    525 CoinBigIndex 
    526 ClpNetworkMatrix::countBasis( const int * whichColumn, 
    527                              int & numberColumnBasic)
    528 {
    529   int i;
    530   CoinBigIndex numberElements=0;
    531   if (trueNetwork_) {
    532     numberElements = 2*numberColumnBasic;
    533   } else {
    534     for (i=0;i<numberColumnBasic;i++) {
    535       int iColumn = whichColumn[i];
    536       CoinBigIndex j=iColumn<<1;
    537       int iRowM = indices_[j];
    538       int iRowP = indices_[j+1];
    539       if (iRowM>=0)
    540         numberElements ++;
    541       if (iRowP>=0)
    542         numberElements ++;
    543     }
    544   }
    545   return numberElements;
     525CoinBigIndex
     526ClpNetworkMatrix::countBasis( const int * whichColumn,
     527                              int & numberColumnBasic)
     528{
     529     int i;
     530     CoinBigIndex numberElements = 0;
     531     if (trueNetwork_) {
     532          numberElements = 2 * numberColumnBasic;
     533     } else {
     534          for (i = 0; i < numberColumnBasic; i++) {
     535               int iColumn = whichColumn[i];
     536               CoinBigIndex j = iColumn << 1;
     537               int iRowM = indices_[j];
     538               int iRowP = indices_[j+1];
     539               if (iRowM >= 0)
     540                    numberElements ++;
     541               if (iRowP >= 0)
     542                    numberElements ++;
     543          }
     544     }
     545     return numberElements;
    546546}
    547547void
    548548ClpNetworkMatrix::fillBasis(ClpSimplex * /*model*/,
    549                          const int * whichColumn,
    550                         int & numberColumnBasic,
    551                         int * indexRowU, int * start,
    552                         int * rowCount, int * columnCount,
    553                         CoinFactorizationDouble * elementU)
    554 {
    555   int i;
    556   CoinBigIndex numberElements=start[0];
    557   if (trueNetwork_) {
    558     for (i=0;i<numberColumnBasic;i++) {
    559       int iColumn = whichColumn[i];
    560       CoinBigIndex j=iColumn<<1;
    561       int iRowM = indices_[j];
    562       int iRowP = indices_[j+1];
    563       indexRowU[numberElements]=iRowM;
    564       rowCount[iRowM]++;
    565       elementU[numberElements]=-1.0;
    566       indexRowU[numberElements+1]=iRowP;
    567       rowCount[iRowP]++;
    568       elementU[numberElements+1]=1.0;
    569       numberElements+=2;
    570       start[i+1]=numberElements;
    571       columnCount[i]=2;
    572     }
    573   } else {
    574     for (i=0;i<numberColumnBasic;i++) {
    575       int iColumn = whichColumn[i];
    576       CoinBigIndex j=iColumn<<1;
    577       int iRowM = indices_[j];
    578       int iRowP = indices_[j+1];
    579       if (iRowM>=0) {
    580         indexRowU[numberElements]=iRowM;
    581         rowCount[iRowM]++;
    582         elementU[numberElements++]=-1.0;
    583       }
    584       if (iRowP>=0) {
    585         indexRowU[numberElements]=iRowP;
    586         rowCount[iRowP]++;
    587         elementU[numberElements++]=1.0;
    588       }
    589       start[i+1]=numberElements;
    590       columnCount[i]=numberElements-start[i];
    591     }
    592   }
     549                            const int * whichColumn,
     550                            int & numberColumnBasic,
     551                            int * indexRowU, int * start,
     552                            int * rowCount, int * columnCount,
     553                            CoinFactorizationDouble * elementU)
     554{
     555     int i;
     556     CoinBigIndex numberElements = start[0];
     557     if (trueNetwork_) {
     558          for (i = 0; i < numberColumnBasic; i++) {
     559               int iColumn = whichColumn[i];
     560               CoinBigIndex j = iColumn << 1;
     561               int iRowM = indices_[j];
     562               int iRowP = indices_[j+1];
     563               indexRowU[numberElements] = iRowM;
     564               rowCount[iRowM]++;
     565               elementU[numberElements] = -1.0;
     566               indexRowU[numberElements+1] = iRowP;
     567               rowCount[iRowP]++;
     568               elementU[numberElements+1] = 1.0;
     569               numberElements += 2;
     570               start[i+1] = numberElements;
     571               columnCount[i] = 2;
     572          }
     573     } else {
     574          for (i = 0; i < numberColumnBasic; i++) {
     575               int iColumn = whichColumn[i];
     576               CoinBigIndex j = iColumn << 1;
     577               int iRowM = indices_[j];
     578               int iRowP = indices_[j+1];
     579               if (iRowM >= 0) {
     580                    indexRowU[numberElements] = iRowM;
     581                    rowCount[iRowM]++;
     582                    elementU[numberElements++] = -1.0;
     583               }
     584               if (iRowP >= 0) {
     585                    indexRowU[numberElements] = iRowP;
     586                    rowCount[iRowP]++;
     587                    elementU[numberElements++] = 1.0;
     588               }
     589               start[i+1] = numberElements;
     590               columnCount[i] = numberElements - start[i];
     591          }
     592     }
    593593}
    594594/* Unpacks a column into an CoinIndexedvector
    595595 */
    596 void 
    597 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);
     596void
     597ClpNetworkMatrix::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);
    607607}
    608608/* Unpacks a column into an CoinIndexedvector
     
    610610Note that model is NOT const.  Bounds and objective could
    611611be modified if doing column generation (just for this variable) */
    612 void 
     612void
    613613ClpNetworkMatrix::unpackPacked(ClpSimplex * /*model*/,
    614                             CoinIndexedVector * rowArray,
    615                             int iColumn) const
    616 {
    617   int * index = rowArray->getIndices();
    618   double * array = rowArray->denseVector();
    619   int number = 0;
    620   CoinBigIndex j=iColumn<<1;
    621   int iRowM = indices_[j];
    622   int iRowP = indices_[j+1];
    623   if (iRowM>=0) {
    624     array[number]=-1.0;
    625     index[number++]=iRowM;
    626   }
    627   if (iRowP>=0) {
    628     array[number]=1.0;
    629     index[number++]=iRowP;
    630   }
    631   rowArray->setNumElements(number);
    632   rowArray->setPackedMode(true);
     614                               CoinIndexedVector * rowArray,
     615                               int iColumn) const
     616{
     617     int * index = rowArray->getIndices();
     618     double * array = rowArray->denseVector();
     619     int number = 0;
     620     CoinBigIndex j = iColumn << 1;
     621     int iRowM = indices_[j];
     622     int iRowP = indices_[j+1];
     623     if (iRowM >= 0) {
     624          array[number] = -1.0;
     625          index[number++] = iRowM;
     626     }
     627     if (iRowP >= 0) {
     628          array[number] = 1.0;
     629          index[number++] = iRowP;
     630     }
     631     rowArray->setNumElements(number);
     632     rowArray->setPackedMode(true);
    633633}
    634634/* Adds multiple of a column into an CoinIndexedvector
    635635      You can use quickAdd to add to vector */
    636 void 
    637 ClpNetworkMatrix::add(const ClpSimplex * /*model*/,CoinIndexedVector * rowArray,
    638                    int iColumn, double multiplier) const
    639 {
    640   CoinBigIndex j=iColumn<<1;
    641   int iRowM = indices_[j];
    642   int iRowP = indices_[j+1];
    643   if (iRowM>=0)
    644     rowArray->quickAdd(iRowM,-multiplier);
    645   if (iRowP>=0)
    646     rowArray->quickAdd(iRowP,multiplier);
     636void
     637ClpNetworkMatrix::add(const ClpSimplex * /*model*/, CoinIndexedVector * rowArray,
     638                      int iColumn, double multiplier) const
     639{
     640     CoinBigIndex j = iColumn << 1;
     641     int iRowM = indices_[j];
     642     int iRowP = indices_[j+1];
     643     if (iRowM >= 0)
     644          rowArray->quickAdd(iRowM, -multiplier);
     645     if (iRowP >= 0)
     646          rowArray->quickAdd(iRowP, multiplier);
    647647}
    648648/* Adds multiple of a column into an array */
    649 void 
    650 ClpNetworkMatrix::add(const ClpSimplex * /*model*/,double * array,
    651                     int iColumn, double multiplier) const
    652 {
    653   CoinBigIndex j=iColumn<<1;
    654   int iRowM = indices_[j];
    655   int iRowP = indices_[j+1];
    656   if (iRowM>=0)
    657     array[iRowM] -= multiplier;
    658   if (iRowP>=0)
    659     array[iRowP] += multiplier;
     649void
     650ClpNetworkMatrix::add(const ClpSimplex * /*model*/, double * array,
     651                      int iColumn, double multiplier) const
     652{
     653     CoinBigIndex j = iColumn << 1;
     654     int iRowM = indices_[j];
     655     int iRowP = indices_[j+1];
     656     if (iRowM >= 0)
     657          array[iRowM] -= multiplier;
     658     if (iRowP >= 0)
     659          array[iRowP] += multiplier;
    660660}
    661661
    662662// Return a complete CoinPackedMatrix
    663 CoinPackedMatrix * 
    664 ClpNetworkMatrix::getPackedMatrix() const 
    665 {
    666   if (!matrix_) {
    667     assert (trueNetwork_); // fix later
    668     int numberElements = 2*numberColumns_;
    669     double * elements = new double [numberElements];
    670     CoinBigIndex i;
    671     for (i=0;i<2*numberColumns_;i+=2) {
    672       elements[i]=-1.0;
    673       elements[i+1]=1.0;
    674     }
    675     CoinBigIndex * starts = new CoinBigIndex [numberColumns_+1];
    676     for (i=0;i<numberColumns_+1;i++) {
    677       starts[i]=2*i;
    678     }
    679     // use assignMatrix to save space
    680     delete [] lengths_;
    681     lengths_ = NULL;
    682     matrix_ = new CoinPackedMatrix();
    683     int * indices = CoinCopyOfArray(indices_,2*numberColumns_);
    684     matrix_->assignMatrix(true,numberRows_,numberColumns_,
    685                           getNumElements(),
    686                           elements,indices,
    687                           starts,lengths_);
    688     assert(!elements);
    689     assert(!starts);
    690     assert (!indices);
    691     assert (!lengths_);
    692   }
    693   return matrix_;
     663CoinPackedMatrix *
     664ClpNetworkMatrix::getPackedMatrix() const
     665{
     666     if (!matrix_) {
     667          assert (trueNetwork_); // fix later
     668          int numberElements = 2 * numberColumns_;
     669          double * elements = new double [numberElements];
     670          CoinBigIndex i;
     671          for (i = 0; i < 2 * numberColumns_; i += 2) {
     672               elements[i] = -1.0;
     673               elements[i+1] = 1.0;
     674          }
     675          CoinBigIndex * starts = new CoinBigIndex [numberColumns_+1];
     676          for (i = 0; i < numberColumns_ + 1; i++) {
     677               starts[i] = 2 * i;
     678          }
     679          // use assignMatrix to save space
     680          delete [] lengths_;
     681          lengths_ = NULL;
     682          matrix_ = new CoinPackedMatrix();
     683          int * indices = CoinCopyOfArray(indices_, 2 * numberColumns_);
     684          matrix_->assignMatrix(true, numberRows_, numberColumns_,
     685                                getNumElements(),
     686                                elements, indices,
     687                                starts, lengths_);
     688          assert(!elements);
     689          assert(!starts);
     690          assert (!indices);
     691          assert (!lengths_);
     692     }
     693     return matrix_;
    694694}
    695695/* A vector containing the elements in the packed matrix. Note that there
     
    697697   major-dimension vector. To get the actual elements one should look at
    698698   this vector together with vectorStarts and vectorLengths. */
    699 const double * 
    700 ClpNetworkMatrix::getElements() const 
    701 {
    702   if (!matrix_)
    703     getPackedMatrix();
    704   return matrix_->getElements();
    705 }
    706 
    707 const CoinBigIndex * 
    708 ClpNetworkMatrix::getVectorStarts() const 
    709 {
    710   if (!matrix_)
    711     getPackedMatrix();
    712   return matrix_->getVectorStarts();
     699const double *
     700ClpNetworkMatrix::getElements() const
     701{
     702     if (!matrix_)
     703          getPackedMatrix();
     704     return matrix_->getElements();
     705}
     706
     707const CoinBigIndex *
     708ClpNetworkMatrix::getVectorStarts() const
     709{
     710     if (!matrix_)
     711          getPackedMatrix();
     712     return matrix_->getVectorStarts();
    713713}
    714714/* The lengths of the major-dimension vectors. */
    715 const int * 
     715const int *
    716716ClpNetworkMatrix::getVectorLengths() const
    717717{
    718   assert (trueNetwork_); // fix later
    719   if (!lengths_) {
    720     lengths_ = new int [numberColumns_];
    721     int i;
    722     for (i=0;i<numberColumns_;i++) {
    723       lengths_[i]=2;
    724     }
    725   }
    726   return lengths_;
     718     assert (trueNetwork_); // fix later
     719     if (!lengths_) {
     720          lengths_ = new int [numberColumns_];
     721          int i;
     722          for (i = 0; i < numberColumns_; i++) {
     723               lengths_[i] = 2;
     724          }
     725     }
     726     return lengths_;
    727727}
    728728/* Delete the columns whose indices are listed in <code>indDel</code>. */
    729 void ClpNetworkMatrix::deleteCols(const int numDel, const int * indDel) 
    730 {
    731   assert (trueNetwork_);
    732   int iColumn;
    733   int numberBad=0;
    734   // Use array to make sure we can have duplicates
    735   char * which = new char[numberColumns_];
    736   memset(which,0,numberColumns_);
    737   int nDuplicate=0;
    738   for (iColumn=0;iColumn<numDel;iColumn++) {
    739     int jColumn = indDel[iColumn];
    740     if (jColumn<0||jColumn>=numberColumns_) {
    741       numberBad++;
    742     } else {
    743       if (which[jColumn])
    744         nDuplicate++;
    745       else
    746         which[jColumn]=1;
    747     }
    748   }
    749   if (numberBad)
    750     throw CoinError("Indices out of range", "deleteCols", "ClpNetworkMatrix");
    751   int newNumber = numberColumns_-numDel+nDuplicate;
    752   // Get rid of temporary arrays
    753   delete [] lengths_;
    754   lengths_=NULL;
    755   delete matrix_;
    756   matrix_= NULL;
    757   int newSize=2*newNumber;
    758   int * newIndices = new int [newSize];
    759   newSize=0;
    760   for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    761     if (!which[iColumn]) {
    762       CoinBigIndex start;
    763       CoinBigIndex i;
    764       start = 2*iColumn;
    765       for (i=start;i<start+2;i++)
    766         newIndices[newSize++]=indices_[i];
    767     }
    768   }
    769   delete [] which;
    770   delete [] indices_;
    771   indices_= newIndices;
    772   numberColumns_ = newNumber;
     729void ClpNetworkMatrix::deleteCols(const int numDel, const int * indDel)
     730{
     731     assert (trueNetwork_);
     732     int iColumn;
     733     int numberBad = 0;
     734     // Use array to make sure we can have duplicates
     735     char * which = new char[numberColumns_];
     736     memset(which, 0, numberColumns_);
     737     int nDuplicate = 0;
     738     for (iColumn = 0; iColumn < numDel; iColumn++) {
     739          int jColumn = indDel[iColumn];
     740          if (jColumn < 0 || jColumn >= numberColumns_) {
     741               numberBad++;
     742          } else {
     743               if (which[jColumn])
     744                    nDuplicate++;
     745               else
     746                    which[jColumn] = 1;
     747          }
     748     }
     749     if (numberBad)
     750          throw CoinError("Indices out of range", "deleteCols", "ClpNetworkMatrix");
     751     int newNumber = numberColumns_ - numDel + nDuplicate;
     752     // Get rid of temporary arrays
     753     delete [] lengths_;
     754     lengths_ = NULL;
     755     delete matrix_;
     756     matrix_ = NULL;
     757     int newSize = 2 * newNumber;
     758     int * newIndices = new int [newSize];
     759     newSize = 0;
     760     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     761          if (!which[iColumn]) {
     762               CoinBigIndex start;
     763               CoinBigIndex i;
     764               start = 2 * iColumn;
     765               for (i = start; i < start + 2; i++)
     766                    newIndices[newSize++] = indices_[i];
     767          }
     768     }
     769     delete [] which;
     770     delete [] indices_;
     771     indices_ = newIndices;
     772     numberColumns_ = newNumber;
    773773}
    774774/* Delete the rows whose indices are listed in <code>indDel</code>. */
    775 void ClpNetworkMatrix::deleteRows(const int numDel, const int * indDel) 
    776 {
    777   int iRow;
    778   int numberBad=0;
    779   // Use array to make sure we can have duplicates
    780   int * which = new int [numberRows_];
    781   memset(which,0,numberRows_*sizeof(int));
    782   for (iRow=0;iRow<numDel;iRow++) {
    783     int jRow = indDel[iRow];
    784     if (jRow<0||jRow>=numberRows_) {
    785       numberBad++;
    786     } else {
    787       which[jRow]=1;
    788     }
    789   }
    790   if (numberBad)
    791     throw CoinError("Indices out of range", "deleteRows", "ClpNetworkMatrix");
    792   // Only valid of all columns have 0 entries
    793   int iColumn;
    794   for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    795     CoinBigIndex start;
    796     CoinBigIndex i;
    797     start = 2*iColumn;
    798     for (i=start;i<start+2;i++) {
    799       int iRow = indices_[i];
    800       if (which[iRow])
    801         numberBad++;
    802     }
    803   }
    804   if (numberBad)
    805     throw CoinError("Row has entries", "deleteRows", "ClpNetworkMatrix");
    806   int newNumber=0;
    807   for (iRow=0;iRow<numberRows_;iRow++) {
    808     if (!which[iRow])
    809       which[iRow]=newNumber++;
    810     else
    811       which[iRow]=-1;
    812   }
    813   for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    814     CoinBigIndex start;
    815     CoinBigIndex i;
    816     start = 2*iColumn;
    817     for (i=start;i<start+2;i++) {
    818       int iRow = indices_[i];
    819       indices_[i]=which[iRow];
    820     }
    821   }
    822   delete [] which;
    823   numberRows_ = newNumber;
     775void ClpNetworkMatrix::deleteRows(const int numDel, const int * indDel)
     776{
     777     int iRow;
     778     int numberBad = 0;
     779     // Use array to make sure we can have duplicates
     780     int * which = new int [numberRows_];
     781     memset(which, 0, numberRows_ * sizeof(int));
     782     for (iRow = 0; iRow < numDel; iRow++) {
     783          int jRow = indDel[iRow];
     784          if (jRow < 0 || jRow >= numberRows_) {
     785               numberBad++;
     786          } else {
     787               which[jRow] = 1;
     788          }
     789     }
     790     if (numberBad)
     791          throw CoinError("Indices out of range", "deleteRows", "ClpNetworkMatrix");
     792     // Only valid of all columns have 0 entries
     793     int iColumn;
     794     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     795          CoinBigIndex start;
     796          CoinBigIndex i;
     797          start = 2 * iColumn;
     798          for (i = start; i < start + 2; i++) {
     799               int iRow = indices_[i];
     800               if (which[iRow])
     801                    numberBad++;
     802          }
     803     }
     804     if (numberBad)
     805          throw CoinError("Row has entries", "deleteRows", "ClpNetworkMatrix");
     806     int newNumber = 0;
     807     for (iRow = 0; iRow < numberRows_; iRow++) {
     808          if (!which[iRow])
     809               which[iRow] = newNumber++;
     810          else
     811               which[iRow] = -1;
     812     }
     813     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     814          CoinBigIndex start;
     815          CoinBigIndex i;
     816          start = 2 * iColumn;
     817          for (i = start; i < start + 2; i++) {
     818               int iRow = indices_[i];
     819               indices_[i] = which[iRow];
     820          }
     821     }
     822     delete [] which;
     823     numberRows_ = newNumber;
    824824}
    825825/* Given positive integer weights for each row fills in sum of weights
     
    827827   Returns weights vector
    828828*/
    829 CoinBigIndex * 
    830 ClpNetworkMatrix::dubiousWeights(const ClpSimplex * model,int * inputWeights) const
    831 {
    832   int numberRows = model->numberRows();
    833   int numberColumns =model->numberColumns();
    834   int number = numberRows+numberColumns;
    835   CoinBigIndex * weights = new CoinBigIndex[number];
    836   int i;
    837   for (i=0;i<numberColumns;i++) {
    838     CoinBigIndex j=i<<1;
    839     CoinBigIndex count=0;
    840     int iRowM = indices_[j];
    841     int iRowP = indices_[j+1];
    842     if (iRowM>=0) {
    843       count += inputWeights[iRowM];
    844     }
    845     if (iRowP>=0) {
    846       count += inputWeights[iRowP];
    847     }
    848     weights[i]=count;
    849   }
    850   for (i=0;i<numberRows;i++) {
    851     weights[i+numberColumns]=inputWeights[i];
    852   }
    853   return weights;
     829CoinBigIndex *
     830ClpNetworkMatrix::dubiousWeights(const ClpSimplex * model, int * inputWeights) const
     831{
     832     int numberRows = model->numberRows();
     833     int numberColumns = model->numberColumns();
     834     int number = numberRows + numberColumns;
     835     CoinBigIndex * weights = new CoinBigIndex[number];
     836     int i;
     837     for (i = 0; i < numberColumns; i++) {
     838          CoinBigIndex j = i << 1;
     839          CoinBigIndex count = 0;
     840          int iRowM = indices_[j];
     841          int iRowP = indices_[j+1];
     842          if (iRowM >= 0) {
     843               count += inputWeights[iRowM];
     844          }
     845          if (iRowP >= 0) {
     846               count += inputWeights[iRowP];
     847          }
     848          weights[i] = count;
     849     }
     850     for (i = 0; i < numberRows; i++) {
     851          weights[i+numberColumns] = inputWeights[i];
     852     }
     853     return weights;
    854854}
    855855/* Returns largest and smallest elements of both signs.
    856856   Largest refers to largest absolute value.
    857857*/
    858 void 
     858void
    859859ClpNetworkMatrix::rangeOfElements(double & smallestNegative, double & largestNegative,
    860                        double & smallestPositive, double & largestPositive)
    861 {
    862   smallestNegative=-1.0;
    863   largestNegative=-1.0;
    864   smallestPositive=1.0;
    865   largestPositive=1.0;
     860                                  double & smallestPositive, double & largestPositive)
     861{
     862     smallestNegative = -1.0;
     863     largestNegative = -1.0;
     864     smallestPositive = 1.0;
     865     largestPositive = 1.0;
    866866}
    867867// Says whether it can do partial pricing
    868 bool 
     868bool
    869869ClpNetworkMatrix::canDoPartialPricing() const
    870870{
    871   return true;
    872 }
    873 // Partial pricing 
    874 void 
     871     return true;
     872}
     873// Partial pricing
     874void
    875875ClpNetworkMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
    876                               int & bestSequence, int & numberWanted)
    877 {
    878   numberWanted=currentWanted_;
    879   int j;
    880   int start = static_cast<int> (startFraction*numberColumns_);
    881   int end = CoinMin(static_cast<int> (endFraction*numberColumns_+1),numberColumns_);
    882   double tolerance=model->currentDualTolerance();
    883   double * reducedCost = model->djRegion();
    884   const double * duals = model->dualRowSolution();
    885   const double * cost = model->costRegion();
    886   double bestDj;
    887   if (bestSequence>=0)
    888     bestDj = fabs(reducedCost[bestSequence]);
    889   else
    890     bestDj=tolerance;
    891   int sequenceOut = model->sequenceOut();
    892   int saveSequence = bestSequence;
    893   if (!trueNetwork_) {
    894     // Not true network
    895     int iSequence;
    896     for (iSequence=start;iSequence<end;iSequence++) {
    897       if (iSequence!=sequenceOut) {
    898         double value;
    899         int iRowM,iRowP;
    900         ClpSimplex::Status status = model->getStatus(iSequence);
    901 
    902         switch(status) {
    903          
    904         case ClpSimplex::basic:
    905         case ClpSimplex::isFixed:
    906           break;
    907         case ClpSimplex::isFree:
    908         case ClpSimplex::superBasic:
    909           value=cost[iSequence];
    910           j = iSequence<<1;
    911           // skip negative rows
    912           iRowM = indices_[j];
    913           iRowP = indices_[j+1];
    914           if (iRowM>=0)
    915             value += duals[iRowM];
    916           if (iRowP>=0)
    917             value -= duals[iRowP];
    918           value = fabs(value);
    919           if (value>FREE_ACCEPT*tolerance) {
    920             numberWanted--;
    921             // we are going to bias towards free (but only if reasonable)
    922             value *= FREE_BIAS;
    923             if (value>bestDj) {
    924               // check flagged variable and correct dj
    925               if (!model->flagged(iSequence)) {
    926                 bestDj=value;
    927                 bestSequence = iSequence;
    928               } else {
    929                 // just to make sure we don't exit before got something
    930                 numberWanted++;
    931               }
    932             }
    933           }
    934           break;
    935         case ClpSimplex::atUpperBound:
    936           value=cost[iSequence];
    937           j = iSequence<<1;
    938           // skip negative rows
    939           iRowM = indices_[j];
    940           iRowP = indices_[j+1];
    941           if (iRowM>=0)
    942             value += duals[iRowM];
    943           if (iRowP>=0)
    944             value -= duals[iRowP];
    945           if (value>tolerance) {
    946             numberWanted--;
    947             if (value>bestDj) {
    948               // check flagged variable and correct dj
    949               if (!model->flagged(iSequence)) {
    950                 bestDj=value;
    951                 bestSequence = iSequence;
    952               } else {
    953                 // just to make sure we don't exit before got something
    954                 numberWanted++;
    955               }
    956             }
    957           }
    958           break;
    959         case ClpSimplex::atLowerBound:
    960           value=cost[iSequence];
    961           j = iSequence<<1;
    962           // skip negative rows
    963           iRowM = indices_[j];
    964           iRowP = indices_[j+1];
    965           if (iRowM>=0)
    966             value += duals[iRowM];
    967           if (iRowP>=0)
    968             value -= duals[iRowP];
    969           value = -value;
    970           if (value>tolerance) {
    971             numberWanted--;
    972             if (value>bestDj) {
    973               // check flagged variable and correct dj
    974               if (!model->flagged(iSequence)) {
    975                 bestDj=value;
    976                 bestSequence = iSequence;
    977               } else {
    978                 // just to make sure we don't exit before got something
    979                 numberWanted++;
    980               }
    981             }
    982           }
    983           break;
    984         }
    985       }
    986       if (!numberWanted)
    987         break;
    988     }
    989     if (bestSequence!=saveSequence) {
    990       // recompute dj
    991       double value=cost[bestSequence];
    992       j = bestSequence<<1;
    993       // skip negative rows
    994       int iRowM = indices_[j];
    995       int iRowP = indices_[j+1];
    996       if (iRowM>=0)
    997         value += duals[iRowM];
    998       if (iRowP>=0)
    999         value -= duals[iRowP];
    1000       reducedCost[bestSequence] = value;
    1001       savedBestSequence_ = bestSequence;
    1002       savedBestDj_ = reducedCost[savedBestSequence_];
    1003     }
    1004   } else {
    1005     // true network
    1006     int iSequence;
    1007     for (iSequence=start;iSequence<end;iSequence++) {
    1008       if (iSequence!=sequenceOut) {
    1009         double value;
    1010         int iRowM,iRowP;
    1011         ClpSimplex::Status status = model->getStatus(iSequence);
    1012 
    1013         switch(status) {
    1014          
    1015         case ClpSimplex::basic:
    1016         case ClpSimplex::isFixed:
    1017           break;
    1018         case ClpSimplex::isFree:
    1019         case ClpSimplex::superBasic:
    1020           value=cost[iSequence];
    1021           j = iSequence<<1;
    1022           iRowM = indices_[j];
    1023           iRowP = indices_[j+1];
    1024           value += duals[iRowM];
    1025           value -= duals[iRowP];
    1026           value = fabs(value);
    1027           if (value>FREE_ACCEPT*tolerance) {
    1028             numberWanted--;
    1029             // we are going to bias towards free (but only if reasonable)
    1030             value *= FREE_BIAS;
    1031             if (value>bestDj) {
    1032               // check flagged variable and correct dj
    1033               if (!model->flagged(iSequence)) {
    1034                 bestDj=value;
    1035                 bestSequence = iSequence;
    1036               } else {
    1037                 // just to make sure we don't exit before got something
    1038                 numberWanted++;
    1039               }
    1040             }
    1041           }
    1042           break;
    1043         case ClpSimplex::atUpperBound:
    1044           value=cost[iSequence];
    1045           j = iSequence<<1;
    1046           iRowM = indices_[j];
    1047           iRowP = indices_[j+1];
    1048           value += duals[iRowM];
    1049           value -= duals[iRowP];
    1050           if (value>tolerance) {
    1051             numberWanted--;
    1052             if (value>bestDj) {
    1053               // check flagged variable and correct dj
    1054               if (!model->flagged(iSequence)) {
    1055                 bestDj=value;
    1056                 bestSequence = iSequence;
    1057               } else {
    1058                 // just to make sure we don't exit before got something
    1059                 numberWanted++;
    1060               }
    1061             }
    1062           }
    1063           break;
    1064         case ClpSimplex::atLowerBound:
    1065           value=cost[iSequence];
    1066           j = iSequence<<1;
    1067           iRowM = indices_[j];
    1068           iRowP = indices_[j+1];
    1069           value += duals[iRowM];
    1070           value -= duals[iRowP];
    1071           value = -value;
    1072           if (value>tolerance) {
    1073             numberWanted--;
    1074             if (value>bestDj) {
    1075               // check flagged variable and correct dj
    1076               if (!model->flagged(iSequence)) {
    1077                 bestDj=value;
    1078                 bestSequence = iSequence;
    1079               } else {
    1080                 // just to make sure we don't exit before got something
    1081                 numberWanted++;
    1082               }
    1083             }
    1084           }
    1085           break;
    1086         }
    1087       }
    1088       if (!numberWanted)
    1089         break;
    1090     }
    1091     if (bestSequence!=saveSequence) {
    1092       // recompute dj
    1093       double value=cost[bestSequence];
    1094       j = bestSequence<<1;
    1095       int iRowM = indices_[j];
    1096       int iRowP = indices_[j+1];
    1097       value += duals[iRowM];
    1098       value -= duals[iRowP];
    1099       reducedCost[bestSequence] = value;
    1100       savedBestSequence_ = bestSequence;
    1101       savedBestDj_ = reducedCost[savedBestSequence_];
    1102     }
    1103   }
    1104   currentWanted_=numberWanted;
     876                                 int & bestSequence, int & numberWanted)
     877{
     878     numberWanted = currentWanted_;
     879     int j;
     880     int start = static_cast<int> (startFraction * numberColumns_);
     881     int end = CoinMin(static_cast<int> (endFraction * numberColumns_ + 1), numberColumns_);
     882     double tolerance = model->currentDualTolerance();
     883     double * reducedCost = model->djRegion();
     884     const double * duals = model->dualRowSolution();
     885     const double * cost = model->costRegion();
     886     double bestDj;
     887     if (bestSequence >= 0)
     888          bestDj = fabs(reducedCost[bestSequence]);
     889     else
     890          bestDj = tolerance;
     891     int sequenceOut = model->sequenceOut();
     892     int saveSequence = bestSequence;
     893     if (!trueNetwork_) {
     894          // Not true network
     895          int iSequence;
     896          for (iSequence = start; iSequence < end; iSequence++) {
     897               if (iSequence != sequenceOut) {
     898                    double value;
     899                    int iRowM, iRowP;
     900                    ClpSimplex::Status status = model->getStatus(iSequence);
     901
     902                    switch(status) {
     903
     904                    case ClpSimplex::basic:
     905                    case ClpSimplex::isFixed:
     906                         break;
     907                    case ClpSimplex::isFree:
     908                    case ClpSimplex::superBasic:
     909                         value = cost[iSequence];
     910                         j = iSequence << 1;
     911                         // skip negative rows
     912                         iRowM = indices_[j];
     913                         iRowP = indices_[j+1];
     914                         if (iRowM >= 0)
     915                              value += duals[iRowM];
     916                         if (iRowP >= 0)
     917                              value -= duals[iRowP];
     918                         value = fabs(value);
     919                         if (value > FREE_ACCEPT * tolerance) {
     920                              numberWanted--;
     921                              // we are going to bias towards free (but only if reasonable)
     922                              value *= FREE_BIAS;
     923                              if (value > bestDj) {
     924                                   // check flagged variable and correct dj
     925                                   if (!model->flagged(iSequence)) {
     926                                        bestDj = value;
     927                                        bestSequence = iSequence;
     928                                   } else {
     929                                        // just to make sure we don't exit before got something
     930                                        numberWanted++;
     931                                   }
     932                              }
     933                         }
     934                         break;
     935                    case ClpSimplex::atUpperBound:
     936                         value = cost[iSequence];
     937                         j = iSequence << 1;
     938                         // skip negative rows
     939                         iRowM = indices_[j];
     940                         iRowP = indices_[j+1];
     941                         if (iRowM >= 0)
     942                              value += duals[iRowM];
     943                         if (iRowP >= 0)
     944                              value -= duals[iRowP];
     945                         if (value > tolerance) {
     946                              numberWanted--;
     947                              if (value > bestDj) {
     948                                   // check flagged variable and correct dj
     949                                   if (!model->flagged(iSequence)) {
     950                                        bestDj = value;
     951                                        bestSequence = iSequence;
     952                                   } else {
     953                                        // just to make sure we don't exit before got something
     954                                        numberWanted++;
     955                                   }
     956                              }
     957                         }
     958                         break;
     959                    case ClpSimplex::atLowerBound:
     960                         value = cost[iSequence];
     961                         j = iSequence << 1;
     962                         // skip negative rows
     963                         iRowM = indices_[j];
     964                         iRowP = indices_[j+1];
     965                         if (iRowM >= 0)
     966                              value += duals[iRowM];
     967                         if (iRowP >= 0)
     968                              value -= duals[iRowP];
     969                         value = -value;
     970                         if (value > tolerance) {
     971                              numberWanted--;
     972                              if (value > bestDj) {
     973                                   // check flagged variable and correct dj
     974                                   if (!model->flagged(iSequence)) {
     975                                        bestDj = value;
     976                                        bestSequence = iSequence;
     977                                   } else {
     978                                        // just to make sure we don't exit before got something
     979                                        numberWanted++;
     980                                   }
     981                              }
     982                         }
     983                         break;
     984                    }
     985               }
     986               if (!numberWanted)
     987                    break;
     988          }
     989          if (bestSequence != saveSequence) {
     990               // recompute dj
     991               double value = cost[bestSequence];
     992               j = bestSequence << 1;
     993               // skip negative rows
     994               int iRowM = indices_[j];
     995               int iRowP = indices_[j+1];
     996               if (iRowM >= 0)
     997                    value += duals[iRowM];
     998               if (iRowP >= 0)
     999                    value -= duals[iRowP];
     1000               reducedCost[bestSequence] = value;
     1001               savedBestSequence_ = bestSequence;
     1002               savedBestDj_ = reducedCost[savedBestSequence_];
     1003          }
     1004     } else {
     1005          // true network
     1006          int iSequence;
     1007          for (iSequence = start; iSequence < end; iSequence++) {
     1008               if (iSequence != sequenceOut) {
     1009                    double value;
     1010                    int iRowM, iRowP;
     1011                    ClpSimplex::Status status = model->getStatus(iSequence);
     1012
     1013                    switch(status) {
     1014
     1015                    case ClpSimplex::basic:
     1016                    case ClpSimplex::isFixed:
     1017                         break;
     1018                    case ClpSimplex::isFree:
     1019                    case ClpSimplex::superBasic:
     1020                         value = cost[iSequence];
     1021                         j = iSequence << 1;
     1022                         iRowM = indices_[j];
     1023                         iRowP = indices_[j+1];
     1024                         value += duals[iRowM];
     1025                         value -= duals[iRowP];
     1026                         value = fabs(value);
     1027                         if (value > FREE_ACCEPT * tolerance) {
     1028                              numberWanted--;
     1029                              // we are going to bias towards free (but only if reasonable)
     1030                              value *= FREE_BIAS;
     1031                              if (value > bestDj) {
     1032                                   // check flagged variable and correct dj
     1033                                   if (!model->flagged(iSequence)) {
     1034                                        bestDj = value;
     1035                                        bestSequence = iSequence;
     1036                                   } else {
     1037                                        // just to make sure we don't exit before got something
     1038                                        numberWanted++;
     1039                                   }
     1040                              }
     1041                         }
     1042                         break;
     1043                    case ClpSimplex::atUpperBound:
     1044                         value = cost[iSequence];
     1045                         j = iSequence << 1;
     1046                         iRowM = indices_[j];
     1047                         iRowP = indices_[j+1];
     1048                         value += duals[iRowM];
     1049                         value -= duals[iRowP];
     1050                         if (value > tolerance) {
     1051                              numberWanted--;
     1052                              if (value > bestDj) {
     1053                                   // check flagged variable and correct dj
     1054                                   if (!model->flagged(iSequence)) {
     1055                                        bestDj = value;
     1056                                        bestSequence = iSequence;
     1057                                   } else {
     1058                                        // just to make sure we don't exit before got something
     1059                                        numberWanted++;
     1060                                   }
     1061                              }
     1062                         }
     1063                         break;
     1064                    case ClpSimplex::atLowerBound:
     1065                         value = cost[iSequence];
     1066                         j = iSequence << 1;
     1067                         iRowM = indices_[j];
     1068                         iRowP = indices_[j+1];
     1069                         value += duals[iRowM];
     1070                         value -= duals[iRowP];
     1071                         value = -value;
     1072                         if (value > tolerance) {
     1073                              numberWanted--;
     1074                              if (value > bestDj) {
     1075                                   // check flagged variable and correct dj
     1076                                   if (!model->flagged(iSequence)) {
     1077                                        bestDj = value;
     1078                                        bestSequence = iSequence;
     1079                                   } else {
     1080                                        // just to make sure we don't exit before got something
     1081                                        numberWanted++;
     1082                                   }
     1083                              }
     1084                         }
     1085                         break;
     1086                    }
     1087               }
     1088               if (!numberWanted)
     1089                    break;
     1090          }
     1091          if (bestSequence != saveSequence) {
     1092               // recompute dj
     1093               double value = cost[bestSequence];
     1094               j = bestSequence << 1;
     1095               int iRowM = indices_[j];
     1096               int iRowP = indices_[j+1];
     1097               value += duals[iRowM];
     1098               value -= duals[iRowP];
     1099               reducedCost[bestSequence] = value;
     1100               savedBestSequence_ = bestSequence;
     1101               savedBestDj_ = reducedCost[savedBestSequence_];
     1102          }
     1103     }
     1104     currentWanted_ = numberWanted;
    11051105}
    11061106// Allow any parts of a created CoinMatrix to be deleted
    1107 void 
    1108 ClpNetworkMatrix::releasePackedMatrix() const 
    1109 {
    1110   delete matrix_;
    1111   delete [] lengths_;
    1112   matrix_=NULL;
    1113   lengths_=NULL;
     1107void
     1108ClpNetworkMatrix::releasePackedMatrix() const
     1109{
     1110     delete matrix_;
     1111     delete [] lengths_;
     1112     matrix_ = NULL;
     1113     lengths_ = NULL;
    11141114}
    11151115// Append Columns
    1116 void 
     1116void
    11171117ClpNetworkMatrix::appendCols(int number, const CoinPackedVectorBase * const * columns)
    11181118{
    1119   int iColumn;
    1120   int numberBad=0;
    1121   for (iColumn=0;iColumn<number;iColumn++) {
    1122     int n=columns[iColumn]->getNumElements();
    1123     const double * element = columns[iColumn]->getElements();
    1124     if (n!=2)
    1125       numberBad++;
    1126     if (fabs(element[0])!=1.0||fabs(element[1])!=1.0)
    1127         numberBad++;
    1128     else if (element[0]*element[1]!=-1.0)
    1129         numberBad++;
    1130   }
    1131   if (numberBad)
    1132     throw CoinError("Not network", "appendCols", "ClpNetworkMatrix");
    1133   // Get rid of temporary arrays
    1134   delete [] lengths_;
    1135   lengths_=NULL;
    1136   delete matrix_;
    1137   matrix_= NULL;
    1138   CoinBigIndex size = 2*number;
    1139   int * temp2 = new int [numberColumns_*2+size];
    1140   CoinMemcpyN(indices_,numberColumns_*2,temp2);
    1141   delete [] indices_;
    1142   indices_= temp2;
    1143   // now add
    1144   size=2*numberColumns_;
    1145   for (iColumn=0;iColumn<number;iColumn++) {
    1146     const int * row = columns[iColumn]->getIndices();
    1147     const double * element = columns[iColumn]->getElements();
    1148     if (element[0]==-1.0) {
    1149       indices_[size++]=row[0];
    1150       indices_[size++]=row[1];
    1151     } else {
    1152       indices_[size++]=row[1];
    1153       indices_[size++]=row[0];
    1154     }
    1155   }
    1156  
    1157   numberColumns_ += number;
     1119     int iColumn;
     1120     int numberBad = 0;
     1121     for (iColumn = 0; iColumn < number; iColumn++) {
     1122          int n = columns[iColumn]->getNumElements();
     1123          const double * element = columns[iColumn]->getElements();
     1124          if (n != 2)
     1125               numberBad++;
     1126          if (fabs(element[0]) != 1.0 || fabs(element[1]) != 1.0)
     1127               numberBad++;
     1128          else if (element[0]*element[1] != -1.0)
     1129               numberBad++;
     1130     }
     1131     if (numberBad)
     1132          throw CoinError("Not network", "appendCols", "ClpNetworkMatrix");
     1133     // Get rid of temporary arrays
     1134     delete [] lengths_;
     1135     lengths_ = NULL;
     1136     delete matrix_;
     1137     matrix_ = NULL;
     1138     CoinBigIndex size = 2 * number;
     1139     int * temp2 = new int [numberColumns_*2+size];
     1140     CoinMemcpyN(indices_, numberColumns_ * 2, temp2);
     1141     delete [] indices_;
     1142     indices_ = temp2;
     1143     // now add
     1144     size = 2 * numberColumns_;
     1145     for (iColumn = 0; iColumn < number; iColumn++) {
     1146          const int * row = columns[iColumn]->getIndices();
     1147          const double * element = columns[iColumn]->getElements();
     1148          if (element[0] == -1.0) {
     1149               indices_[size++] = row[0];
     1150               indices_[size++] = row[1];
     1151          } else {
     1152               indices_[size++] = row[1];
     1153               indices_[size++] = row[0];
     1154          }
     1155     }
     1156
     1157     numberColumns_ += number;
    11581158}
    11591159// Append Rows
    1160 void 
     1160void
    11611161ClpNetworkMatrix::appendRows(int number, const CoinPackedVectorBase * const * rows)
    11621162{
    1163   // must be zero arrays
    1164   int numberBad=0;
    1165   int iRow;
    1166   for (iRow=0;iRow<number;iRow++) {
    1167     numberBad += rows[iRow]->getNumElements();
    1168   }
    1169   if (numberBad)
    1170     throw CoinError("Not NULL rows", "appendRows", "ClpNetworkMatrix");
    1171   numberRows_ += number;
     1163     // must be zero arrays
     1164     int numberBad = 0;
     1165     int iRow;
     1166     for (iRow = 0; iRow < number; iRow++) {
     1167          numberBad += rows[iRow]->getNumElements();
     1168     }
     1169     if (numberBad)
     1170          throw CoinError("Not NULL rows", "appendRows", "ClpNetworkMatrix");
     1171     numberRows_ += number;
    11721172}
    11731173#ifndef SLIM_CLP
     
    11761176   number of columns-1/rows-1 (if numberOther>0) or duplicates
    11771177   If 0 then rows, 1 if columns */
    1178 int 
     1178int
    11791179ClpNetworkMatrix::appendMatrix(int number, int type,
    1180                                     const CoinBigIndex * starts, const int * index,
    1181                                const double * element, int /*numberOther*/)
    1182 {
    1183   int numberErrors=0;
    1184   // make into CoinPackedVector
    1185   CoinPackedVectorBase ** vectors=
    1186     new CoinPackedVectorBase * [number];
    1187   int iVector;
    1188   for (iVector=0;iVector<number;iVector++) {
    1189     int iStart = starts[iVector];
    1190     vectors[iVector] =
    1191       new CoinPackedVector(starts[iVector+1]-iStart,
    1192                            index+iStart,element+iStart);
    1193   }
    1194   if (type==0) {
    1195     // rows
    1196     appendRows(number,vectors);
    1197   } else {
    1198     // columns
    1199     appendCols(number,vectors);
    1200   }
    1201   for (iVector=0;iVector<number;iVector++)
    1202     delete vectors[iVector];
    1203   delete [] vectors;
    1204   return numberErrors;
     1180                               const CoinBigIndex * starts, const int * index,
     1181                               const double * element, int /*numberOther*/)
     1182{
     1183     int numberErrors = 0;
     1184     // make into CoinPackedVector
     1185     CoinPackedVectorBase ** vectors =
     1186          new CoinPackedVectorBase * [number];
     1187     int iVector;
     1188     for (iVector = 0; iVector < number; iVector++) {
     1189          int iStart = starts[iVector];
     1190          vectors[iVector] =
     1191               new CoinPackedVector(starts[iVector+1] - iStart,
     1192                                    index + iStart, element + iStart);
     1193     }
     1194     if (type == 0) {
     1195          // rows
     1196          appendRows(number, vectors);
     1197     } else {
     1198          // columns
     1199          appendCols(number, vectors);
     1200     }
     1201     for (iVector = 0; iVector < number; iVector++)
     1202          delete vectors[iVector];
     1203     delete [] vectors;
     1204     return numberErrors;
    12051205}
    12061206#endif
    12071207/* Subset clone (without gaps).  Duplicates are allowed
    12081208   and order is as given */
    1209 ClpMatrixBase * 
     1209ClpMatrixBase *
    12101210ClpNetworkMatrix::subsetClone (int numberRows, const int * whichRows,
    1211                               int numberColumns,
    1212                               const int * whichColumns) const
    1213 {
    1214   return new ClpNetworkMatrix(*this, numberRows, whichRows,
    1215                                    numberColumns, whichColumns);
     1211                               int numberColumns,
     1212                               const int * whichColumns) const
     1213{
     1214     return new ClpNetworkMatrix(*this, numberRows, whichRows,
     1215                                 numberColumns, whichColumns);
    12161216}
    12171217/* Subset constructor (without gaps).  Duplicates are allowed
    12181218   and order is as given */
    12191219ClpNetworkMatrix::ClpNetworkMatrix (
    1220                        const ClpNetworkMatrix & rhs,
    1221                        int numberRows, const int * whichRow,
    1222                        int numberColumns, const int * whichColumn)
    1223 : ClpMatrixBase(rhs)
    1224 { 
    1225   setType(11);
    1226   matrix_ = NULL;
    1227   lengths_=NULL;
    1228   indices_=new int[2*numberColumns];;
    1229   numberRows_=numberRows;
    1230   numberColumns_=numberColumns;
    1231   trueNetwork_=true;
    1232   int iColumn;
    1233   int numberBad=0;
    1234   int * which = new int [rhs.numberRows_];
    1235   int iRow;
    1236   for (iRow=0;iRow<rhs.numberRows_;iRow++)
    1237     which[iRow]=-1;
    1238   int n=0;
    1239   for (iRow=0;iRow<numberRows;iRow++) {
    1240     int jRow=whichRow[iRow];
    1241     assert (jRow>=0&&jRow<rhs.numberRows_);
    1242     which[jRow]=n++;
    1243   }
    1244   for (iColumn=0;iColumn<numberColumns;iColumn++) {
    1245     CoinBigIndex start;
    1246     CoinBigIndex i;
    1247     start = 2*iColumn;
    1248     CoinBigIndex offset = 2*whichColumn[iColumn]-start;
    1249     for (i=start;i<start+2;i++) {
    1250       int iRow = rhs.indices_[i+offset];
    1251       iRow = which[iRow];
    1252       if (iRow<0)
    1253         numberBad++;
    1254       else
    1255         indices_[i]=iRow;
    1256     }
    1257   }
    1258   if (numberBad)
    1259     throw CoinError("Invalid rows", "subsetConstructor", "ClpNetworkMatrix");
    1260 }
     1220     const ClpNetworkMatrix & rhs,
     1221     int numberRows, const int * whichRow,
     1222     int numberColumns, const int * whichColumn)
     1223     : ClpMatrixBase(rhs)
     1224{
     1225     setType(11);
     1226     matrix_ = NULL;
     1227     lengths_ = NULL;
     1228     indices_ = new int[2*numberColumns];;
     1229     numberRows_ = numberRows;
     1230     numberColumns_ = numberColumns;
     1231     trueNetwork_ = true;
     1232     int iColumn;
     1233     int numberBad = 0;
     1234     int * which = new int [rhs.numberRows_];
     1235     int iRow;
     1236     for (iRow = 0; iRow < rhs.numberRows_; iRow++)
     1237          which[iRow] = -1;
     1238     int n = 0;
     1239     for (iRow = 0; iRow < numberRows; iRow++) {
     1240          int jRow = whichRow[iRow];
     1241          assert (jRow >= 0 && jRow < rhs.numberRows_);
     1242          which[jRow] = n++;
     1243     }
     1244     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     1245          CoinBigIndex start;
     1246          CoinBigIndex i;
     1247          start = 2 * iColumn;
     1248          CoinBigIndex offset = 2 * whichColumn[iColumn] - start;
     1249          for (i = start; i < start + 2; i++) {
     1250               int iRow = rhs.indices_[i+offset];
     1251               iRow = which[iRow];
     1252               if (iRow < 0)
     1253                    numberBad++;
     1254               else
     1255                    indices_[i] = iRow;
     1256          }
     1257     }
     1258     if (numberBad)
     1259          throw CoinError("Invalid rows", "subsetConstructor", "ClpNetworkMatrix");
     1260}
Note: See TracChangeset for help on using the changeset viewer.