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

formatting

File:
1 edited

Legend:

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

    r2070 r2385  
    1818// Default Constructor
    1919//-------------------------------------------------------------------
    20 ClpDualRowSteepest::ClpDualRowSteepest (int mode)
    21      : ClpDualRowPivot(),
    22        state_(-1),
    23        mode_(mode),
    24        persistence_(normal),
    25        weights_(NULL),
    26        infeasible_(NULL),
    27        alternateWeights_(NULL),
    28        savedWeights_(NULL),
    29       dubiousWeights_(NULL)
    30 {
    31      type_ = 2 + 64 * mode;
     20ClpDualRowSteepest::ClpDualRowSteepest(int mode)
     21  : ClpDualRowPivot()
     22  , state_(-1)
     23  , mode_(mode)
     24  , persistence_(normal)
     25  , weights_(NULL)
     26  , infeasible_(NULL)
     27  , alternateWeights_(NULL)
     28  , savedWeights_(NULL)
     29  , dubiousWeights_(NULL)
     30{
     31  type_ = 2 + 64 * mode;
    3232}
    3333
     
    3535// Copy constructor
    3636//-------------------------------------------------------------------
    37 ClpDualRowSteepest::ClpDualRowSteepest (const ClpDualRowSteepest & rhs)
    38      : ClpDualRowPivot(rhs)
    39 {
    40      state_ = rhs.state_;
    41      mode_ = rhs.mode_;
    42      persistence_ = rhs.persistence_;
    43      model_ = rhs.model_;
    44      if ((model_ && model_->whatsChanged() & 1) != 0) {
    45           int number = model_->numberRows();
    46           if (rhs.savedWeights_)
    47                number = CoinMin(number, rhs.savedWeights_->capacity());
    48           if (rhs.infeasible_) {
    49                infeasible_ = new CoinIndexedVector(rhs.infeasible_);
    50           } else {
    51                infeasible_ = NULL;
    52           }
    53           if (rhs.weights_) {
    54                weights_ = new double[number];
    55                ClpDisjointCopyN(rhs.weights_, number, weights_);
    56           } else {
    57                weights_ = NULL;
    58           }
    59           if (rhs.alternateWeights_) {
    60                alternateWeights_ = new CoinIndexedVector(rhs.alternateWeights_);
    61           } else {
    62                alternateWeights_ = NULL;
    63           }
    64           if (rhs.savedWeights_) {
    65                savedWeights_ = new CoinIndexedVector(rhs.savedWeights_);
    66           } else {
    67                savedWeights_ = NULL;
    68           }
    69           if (rhs.dubiousWeights_) {
    70                assert(model_);
    71                int number = model_->numberRows();
    72                dubiousWeights_ = new int[number];
    73                ClpDisjointCopyN(rhs.dubiousWeights_, number, dubiousWeights_);
    74           } else {
    75                dubiousWeights_ = NULL;
    76           }
    77      } else {
    78           infeasible_ = NULL;
    79           weights_ = NULL;
    80           alternateWeights_ = NULL;
    81           savedWeights_ = NULL;
    82           dubiousWeights_ = NULL;
    83      }
     37ClpDualRowSteepest::ClpDualRowSteepest(const ClpDualRowSteepest &rhs)
     38  : ClpDualRowPivot(rhs)
     39{
     40  state_ = rhs.state_;
     41  mode_ = rhs.mode_;
     42  persistence_ = rhs.persistence_;
     43  model_ = rhs.model_;
     44  if ((model_ && model_->whatsChanged() & 1) != 0) {
     45    int number = model_->numberRows();
     46    if (rhs.savedWeights_)
     47      number = CoinMin(number, rhs.savedWeights_->capacity());
     48    if (rhs.infeasible_) {
     49      infeasible_ = new CoinIndexedVector(rhs.infeasible_);
     50    } else {
     51      infeasible_ = NULL;
     52    }
     53    if (rhs.weights_) {
     54      weights_ = new double[number];
     55      ClpDisjointCopyN(rhs.weights_, number, weights_);
     56    } else {
     57      weights_ = NULL;
     58    }
     59    if (rhs.alternateWeights_) {
     60      alternateWeights_ = new CoinIndexedVector(rhs.alternateWeights_);
     61    } else {
     62      alternateWeights_ = NULL;
     63    }
     64    if (rhs.savedWeights_) {
     65      savedWeights_ = new CoinIndexedVector(rhs.savedWeights_);
     66    } else {
     67      savedWeights_ = NULL;
     68    }
     69    if (rhs.dubiousWeights_) {
     70      assert(model_);
     71      int number = model_->numberRows();
     72      dubiousWeights_ = new int[number];
     73      ClpDisjointCopyN(rhs.dubiousWeights_, number, dubiousWeights_);
     74    } else {
     75      dubiousWeights_ = NULL;
     76    }
     77  } else {
     78    infeasible_ = NULL;
     79    weights_ = NULL;
     80    alternateWeights_ = NULL;
     81    savedWeights_ = NULL;
     82    dubiousWeights_ = NULL;
     83  }
    8484}
    8585
     
    8787// Destructor
    8888//-------------------------------------------------------------------
    89 ClpDualRowSteepest::~ClpDualRowSteepest ()
    90 {
    91      delete [] weights_;
    92      delete [] dubiousWeights_;
    93      delete infeasible_;
    94      delete alternateWeights_;
    95      delete savedWeights_;
     89ClpDualRowSteepest::~ClpDualRowSteepest()
     90{
     91  delete[] weights_;
     92  delete[] dubiousWeights_;
     93  delete infeasible_;
     94  delete alternateWeights_;
     95  delete savedWeights_;
    9696}
    9797
     
    100100//-------------------------------------------------------------------
    101101ClpDualRowSteepest &
    102 ClpDualRowSteepest::operator=(const ClpDualRowSteepest& rhs)
    103 {
    104      if (this != &rhs) {
    105           ClpDualRowPivot::operator=(rhs);
    106           state_ = rhs.state_;
    107           mode_ = rhs.mode_;
    108           persistence_ = rhs.persistence_;
    109           model_ = rhs.model_;
    110           delete [] weights_;
    111           delete [] dubiousWeights_;
    112           delete infeasible_;
    113           delete alternateWeights_;
    114           delete savedWeights_;
    115           assert(model_);
    116           int number = model_->numberRows();
    117           if (rhs.savedWeights_)
    118                number = CoinMin(number, rhs.savedWeights_->capacity());
    119           if (rhs.infeasible_ != NULL) {
    120                infeasible_ = new CoinIndexedVector(rhs.infeasible_);
     102ClpDualRowSteepest::operator=(const ClpDualRowSteepest &rhs)
     103{
     104  if (this != &rhs) {
     105    ClpDualRowPivot::operator=(rhs);
     106    state_ = rhs.state_;
     107    mode_ = rhs.mode_;
     108    persistence_ = rhs.persistence_;
     109    model_ = rhs.model_;
     110    delete[] weights_;
     111    delete[] dubiousWeights_;
     112    delete infeasible_;
     113    delete alternateWeights_;
     114    delete savedWeights_;
     115    assert(model_);
     116    int number = model_->numberRows();
     117    if (rhs.savedWeights_)
     118      number = CoinMin(number, rhs.savedWeights_->capacity());
     119    if (rhs.infeasible_ != NULL) {
     120      infeasible_ = new CoinIndexedVector(rhs.infeasible_);
     121    } else {
     122      infeasible_ = NULL;
     123    }
     124    if (rhs.weights_ != NULL) {
     125      weights_ = new double[number];
     126      ClpDisjointCopyN(rhs.weights_, number, weights_);
     127    } else {
     128      weights_ = NULL;
     129    }
     130    if (rhs.alternateWeights_ != NULL) {
     131      alternateWeights_ = new CoinIndexedVector(rhs.alternateWeights_);
     132    } else {
     133      alternateWeights_ = NULL;
     134    }
     135    if (rhs.savedWeights_ != NULL) {
     136      savedWeights_ = new CoinIndexedVector(rhs.savedWeights_);
     137    } else {
     138      savedWeights_ = NULL;
     139    }
     140    if (rhs.dubiousWeights_) {
     141      assert(model_);
     142      int number = model_->numberRows();
     143      dubiousWeights_ = new int[number];
     144      ClpDisjointCopyN(rhs.dubiousWeights_, number, dubiousWeights_);
     145    } else {
     146      dubiousWeights_ = NULL;
     147    }
     148  }
     149  return *this;
     150}
     151// Fill most values
     152void ClpDualRowSteepest::fill(const ClpDualRowSteepest &rhs)
     153{
     154  state_ = rhs.state_;
     155  mode_ = rhs.mode_;
     156  persistence_ = rhs.persistence_;
     157  assert(model_->numberRows() == rhs.model_->numberRows());
     158  model_ = rhs.model_;
     159  assert(model_);
     160  int number = model_->numberRows();
     161  if (rhs.savedWeights_)
     162    number = CoinMin(number, rhs.savedWeights_->capacity());
     163  if (rhs.infeasible_ != NULL) {
     164    if (!infeasible_)
     165      infeasible_ = new CoinIndexedVector(rhs.infeasible_);
     166    else
     167      *infeasible_ = *rhs.infeasible_;
     168  } else {
     169    delete infeasible_;
     170    infeasible_ = NULL;
     171  }
     172  if (rhs.weights_ != NULL) {
     173    if (!weights_)
     174      weights_ = new double[number];
     175    ClpDisjointCopyN(rhs.weights_, number, weights_);
     176  } else {
     177    delete[] weights_;
     178    weights_ = NULL;
     179  }
     180  if (rhs.alternateWeights_ != NULL) {
     181    if (!alternateWeights_)
     182      alternateWeights_ = new CoinIndexedVector(rhs.alternateWeights_);
     183    else
     184      *alternateWeights_ = *rhs.alternateWeights_;
     185  } else {
     186    delete alternateWeights_;
     187    alternateWeights_ = NULL;
     188  }
     189  if (rhs.savedWeights_ != NULL) {
     190    if (!savedWeights_)
     191      savedWeights_ = new CoinIndexedVector(rhs.savedWeights_);
     192    else
     193      *savedWeights_ = *rhs.savedWeights_;
     194  } else {
     195    delete savedWeights_;
     196    savedWeights_ = NULL;
     197  }
     198  if (rhs.dubiousWeights_) {
     199    assert(model_);
     200    int number = model_->numberRows();
     201    if (!dubiousWeights_)
     202      dubiousWeights_ = new int[number];
     203    ClpDisjointCopyN(rhs.dubiousWeights_, number, dubiousWeights_);
     204  } else {
     205    delete[] dubiousWeights_;
     206    dubiousWeights_ = NULL;
     207  }
     208}
     209// Returns pivot row, -1 if none
     210int ClpDualRowSteepest::pivotRow()
     211{
     212  assert(model_);
     213  int i, iRow;
     214  double *infeas = infeasible_->denseVector();
     215  double largest = 0.0;
     216  int *index = infeasible_->getIndices();
     217  int number = infeasible_->getNumElements();
     218  const int *pivotVariable = model_->pivotVariable();
     219  int chosenRow = -1;
     220  int lastPivotRow = model_->pivotRow();
     221  assert(lastPivotRow < model_->numberRows());
     222  double tolerance = model_->currentPrimalTolerance();
     223  // we can't really trust infeasibilities if there is primal error
     224  // this coding has to mimic coding in checkPrimalSolution
     225  double error = CoinMin(1.0e-2, model_->largestPrimalError());
     226  // allow tolerance at least slightly bigger than standard
     227  tolerance = tolerance + error;
     228  // But cap
     229  tolerance = CoinMin(1000.0, tolerance);
     230  tolerance *= tolerance; // as we are using squares
     231  bool toleranceChanged = false;
     232  double *solution = model_->solutionRegion();
     233  double *lower = model_->lowerRegion();
     234  double *upper = model_->upperRegion();
     235  // do last pivot row one here
     236  //#define CLP_DUAL_FIXED_COLUMN_MULTIPLIER 10.0
     237  if (lastPivotRow >= 0 && lastPivotRow < model_->numberRows()) {
     238#ifdef CLP_DUAL_COLUMN_MULTIPLIER
     239    int numberColumns = model_->numberColumns();
     240#endif
     241    int iPivot = pivotVariable[lastPivotRow];
     242    double value = solution[iPivot];
     243    double lower = model_->lower(iPivot);
     244    double upper = model_->upper(iPivot);
     245    if (value > upper + tolerance) {
     246      value -= upper;
     247      value *= value;
     248#ifdef CLP_DUAL_COLUMN_MULTIPLIER
     249      if (iPivot < numberColumns)
     250        value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
     251#endif
     252      // store square in list
     253      if (infeas[lastPivotRow])
     254        infeas[lastPivotRow] = value; // already there
     255      else
     256        infeasible_->quickAdd(lastPivotRow, value);
     257    } else if (value < lower - tolerance) {
     258      value -= lower;
     259      value *= value;
     260#ifdef CLP_DUAL_COLUMN_MULTIPLIER
     261      if (iPivot < numberColumns)
     262        value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
     263#endif
     264      // store square in list
     265      if (infeas[lastPivotRow])
     266        infeas[lastPivotRow] = value; // already there
     267      else
     268        infeasible_->add(lastPivotRow, value);
     269    } else {
     270      // feasible - was it infeasible - if so set tiny
     271      if (infeas[lastPivotRow])
     272        infeas[lastPivotRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
     273    }
     274    number = infeasible_->getNumElements();
     275  }
     276  if (model_->numberIterations() < model_->lastBadIteration() + 200) {
     277    // we can't really trust infeasibilities if there is dual error
     278    if (model_->largestDualError() > model_->largestPrimalError()) {
     279      tolerance *= CoinMin(model_->largestDualError() / model_->largestPrimalError(), 1000.0);
     280      toleranceChanged = true;
     281    }
     282  }
     283  int numberWanted;
     284  if (mode_ < 2) {
     285    numberWanted = number + 1;
     286  } else if (mode_ == 2) {
     287    numberWanted = CoinMax(2000, number / 8);
     288  } else {
     289    int numberElements = model_->factorization()->numberElements();
     290    double ratio = static_cast< double >(numberElements) / static_cast< double >(model_->numberRows());
     291    numberWanted = CoinMax(2000, number / 8);
     292    if (ratio < 1.0) {
     293      numberWanted = CoinMax(2000, number / 20);
     294    } else if (ratio > 10.0) {
     295      ratio = number * (ratio / 80.0);
     296      if (ratio > number)
     297        numberWanted = number + 1;
     298      else
     299        numberWanted = CoinMax(2000, static_cast< int >(ratio));
     300    }
     301  }
     302  if (model_->largestPrimalError() > 1.0e-3)
     303    numberWanted = number + 1; // be safe
     304  int iPass;
     305  // Setup two passes
     306  int start[4];
     307  start[1] = number;
     308  start[2] = 0;
     309  double dstart = static_cast< double >(number) * model_->randomNumberGenerator()->randomDouble();
     310  start[0] = static_cast< int >(dstart);
     311  start[3] = start[0];
     312  //double largestWeight=0.0;
     313  //double smallestWeight=1.0e100;
     314  for (iPass = 0; iPass < 2; iPass++) {
     315    int end = start[2 * iPass + 1];
     316    for (i = start[2 * iPass]; i < end; i++) {
     317      iRow = index[i];
     318      double value = infeas[iRow];
     319      if (value > tolerance) {
     320        //#define OUT_EQ
     321#ifdef OUT_EQ
     322        {
     323          int iSequence = pivotVariable[iRow];
     324          if (upper[iSequence] == lower[iSequence])
     325            value *= 2.0;
     326        }
     327#endif
     328        double weight = CoinMin(weights_[iRow], 1.0e50);
     329        //largestWeight = CoinMax(largestWeight,weight);
     330        //smallestWeight = CoinMin(smallestWeight,weight);
     331        //double dubious = dubiousWeights_[iRow];
     332        //weight *= dubious;
     333        //if (value>2.0*largest*weight||(value>0.5*largest*weight&&value*largestWeight>dubious*largest*weight)) {
     334        if (value > largest * weight) {
     335          // make last pivot row last resort choice
     336          if (iRow == lastPivotRow) {
     337            if (value * 1.0e-10 < largest * weight)
     338              continue;
     339            else
     340              value *= 1.0e-10;
     341          }
     342          int iSequence = pivotVariable[iRow];
     343          if (!model_->flagged(iSequence)) {
     344            //#define CLP_DEBUG 3
     345#ifdef CLP_DEBUG
     346            double value2 = 0.0;
     347            if (solution[iSequence] > upper[iSequence] + tolerance)
     348              value2 = solution[iSequence] - upper[iSequence];
     349            else if (solution[iSequence] < lower[iSequence] - tolerance)
     350              value2 = solution[iSequence] - lower[iSequence];
     351            assert(fabs(value2 * value2 - infeas[iRow]) < 1.0e-8 * CoinMin(value2 * value2, infeas[iRow]));
     352#endif
     353            if (solution[iSequence] > upper[iSequence] + tolerance || solution[iSequence] < lower[iSequence] - tolerance) {
     354              chosenRow = iRow;
     355              largest = value / weight;
     356              //largestWeight = dubious;
     357            }
    121358          } else {
    122                infeasible_ = NULL;
     359            // just to make sure we don't exit before got something
     360            numberWanted++;
    123361          }
    124           if (rhs.weights_ != NULL) {
    125                weights_ = new double[number];
    126                ClpDisjointCopyN(rhs.weights_, number, weights_);
    127           } else {
    128                weights_ = NULL;
    129           }
    130           if (rhs.alternateWeights_ != NULL) {
    131                alternateWeights_ = new CoinIndexedVector(rhs.alternateWeights_);
    132           } else {
    133                alternateWeights_ = NULL;
    134           }
    135           if (rhs.savedWeights_ != NULL) {
    136                savedWeights_ = new CoinIndexedVector(rhs.savedWeights_);
    137           } else {
    138                savedWeights_ = NULL;
    139           }
    140           if (rhs.dubiousWeights_) {
    141                assert(model_);
    142                int number = model_->numberRows();
    143                dubiousWeights_ = new int[number];
    144                ClpDisjointCopyN(rhs.dubiousWeights_, number, dubiousWeights_);
    145           } else {
    146                dubiousWeights_ = NULL;
    147           }
    148      }
    149      return *this;
    150 }
    151 // Fill most values
    152 void
    153 ClpDualRowSteepest::fill(const ClpDualRowSteepest& rhs)
    154 {
    155      state_ = rhs.state_;
    156      mode_ = rhs.mode_;
    157      persistence_ = rhs.persistence_;
    158      assert (model_->numberRows() == rhs.model_->numberRows());
    159      model_ = rhs.model_;
    160      assert(model_);
    161      int number = model_->numberRows();
    162      if (rhs.savedWeights_)
    163           number = CoinMin(number, rhs.savedWeights_->capacity());
    164      if (rhs.infeasible_ != NULL) {
    165           if (!infeasible_)
    166                infeasible_ = new CoinIndexedVector(rhs.infeasible_);
    167           else
    168                *infeasible_ = *rhs.infeasible_;
    169      } else {
    170           delete infeasible_;
    171           infeasible_ = NULL;
    172      }
    173      if (rhs.weights_ != NULL) {
    174           if (!weights_)
    175                weights_ = new double[number];
    176           ClpDisjointCopyN(rhs.weights_, number, weights_);
    177      } else {
    178           delete [] weights_;
    179           weights_ = NULL;
    180      }
    181      if (rhs.alternateWeights_ != NULL) {
    182           if (!alternateWeights_)
    183                alternateWeights_ = new CoinIndexedVector(rhs.alternateWeights_);
    184           else
    185                *alternateWeights_ = *rhs.alternateWeights_;
    186      } else {
    187           delete alternateWeights_;
    188           alternateWeights_ = NULL;
    189      }
    190      if (rhs.savedWeights_ != NULL) {
    191           if (!savedWeights_)
    192                savedWeights_ = new CoinIndexedVector(rhs.savedWeights_);
    193           else
    194                *savedWeights_ = *rhs.savedWeights_;
    195      } else {
    196           delete savedWeights_;
    197           savedWeights_ = NULL;
    198      }
    199      if (rhs.dubiousWeights_) {
    200           assert(model_);
    201           int number = model_->numberRows();
    202           if (!dubiousWeights_)
    203                dubiousWeights_ = new int[number];
    204           ClpDisjointCopyN(rhs.dubiousWeights_, number, dubiousWeights_);
    205      } else {
    206           delete [] dubiousWeights_;
    207           dubiousWeights_ = NULL;
    208      }
    209 }
    210 // Returns pivot row, -1 if none
    211 int
    212 ClpDualRowSteepest::pivotRow()
    213 {
    214      assert(model_);
    215      int i, iRow;
    216      double * infeas = infeasible_->denseVector();
    217      double largest = 0.0;
    218      int * index = infeasible_->getIndices();
    219      int number = infeasible_->getNumElements();
    220      const int * pivotVariable = model_->pivotVariable();
    221      int chosenRow = -1;
    222      int lastPivotRow = model_->pivotRow();
    223      assert (lastPivotRow < model_->numberRows());
    224      double tolerance = model_->currentPrimalTolerance();
    225      // we can't really trust infeasibilities if there is primal error
    226      // this coding has to mimic coding in checkPrimalSolution
    227      double error = CoinMin(1.0e-2, model_->largestPrimalError());
    228      // allow tolerance at least slightly bigger than standard
    229      tolerance = tolerance +  error;
    230      // But cap
    231      tolerance = CoinMin(1000.0, tolerance);
    232      tolerance *= tolerance; // as we are using squares
    233      bool toleranceChanged = false;
    234      double * solution = model_->solutionRegion();
    235      double * lower = model_->lowerRegion();
    236      double * upper = model_->upperRegion();
    237      // do last pivot row one here
    238      //#define CLP_DUAL_FIXED_COLUMN_MULTIPLIER 10.0
    239      if (lastPivotRow >= 0 && lastPivotRow < model_->numberRows()) {
    240 #ifdef CLP_DUAL_COLUMN_MULTIPLIER
    241           int numberColumns = model_->numberColumns();
    242 #endif
    243           int iPivot = pivotVariable[lastPivotRow];
    244           double value = solution[iPivot];
    245           double lower = model_->lower(iPivot);
    246           double upper = model_->upper(iPivot);
    247           if (value > upper + tolerance) {
    248                value -= upper;
    249                value *= value;
    250 #ifdef CLP_DUAL_COLUMN_MULTIPLIER
    251                if (iPivot < numberColumns)
    252                     value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
    253 #endif
    254                // store square in list
    255                if (infeas[lastPivotRow])
    256                     infeas[lastPivotRow] = value; // already there
    257                else
    258                     infeasible_->quickAdd(lastPivotRow, value);
    259           } else if (value < lower - tolerance) {
    260                value -= lower;
    261                value *= value;
    262 #ifdef CLP_DUAL_COLUMN_MULTIPLIER
    263                if (iPivot < numberColumns)
    264                     value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
    265 #endif
    266                // store square in list
    267                if (infeas[lastPivotRow])
    268                     infeas[lastPivotRow] = value; // already there
    269                else
    270                     infeasible_->add(lastPivotRow, value);
    271           } else {
    272                // feasible - was it infeasible - if so set tiny
    273                if (infeas[lastPivotRow])
    274                     infeas[lastPivotRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
    275           }
    276           number = infeasible_->getNumElements();
    277      }
    278      if(model_->numberIterations() < model_->lastBadIteration() + 200) {
    279           // we can't really trust infeasibilities if there is dual error
    280           if (model_->largestDualError() > model_->largestPrimalError()) {
    281                tolerance *= CoinMin(model_->largestDualError() / model_->largestPrimalError(), 1000.0);
    282                toleranceChanged = true;
    283           }
    284      }
    285      int numberWanted;
    286      if (mode_ < 2 ) {
    287           numberWanted = number + 1;
    288      } else if (mode_ == 2) {
    289           numberWanted = CoinMax(2000, number / 8);
    290      } else {
    291           int numberElements = model_->factorization()->numberElements();
    292           double ratio = static_cast<double> (numberElements) /
    293                          static_cast<double> (model_->numberRows());
    294           numberWanted = CoinMax(2000, number / 8);
    295           if (ratio < 1.0) {
    296                numberWanted = CoinMax(2000, number / 20);
    297           } else if (ratio > 10.0) {
    298                ratio = number * (ratio / 80.0);
    299                if (ratio > number)
    300                     numberWanted = number + 1;
    301                else
    302                     numberWanted = CoinMax(2000, static_cast<int> (ratio));
    303           }
    304      }
    305      if (model_->largestPrimalError() > 1.0e-3)
    306           numberWanted = number + 1; // be safe
    307      int iPass;
    308      // Setup two passes
    309      int start[4];
    310      start[1] = number;
    311      start[2] = 0;
    312      double dstart = static_cast<double> (number) * model_->randomNumberGenerator()->randomDouble();
    313      start[0] = static_cast<int> (dstart);
    314      start[3] = start[0];
    315      //double largestWeight=0.0;
    316      //double smallestWeight=1.0e100;
    317      for (iPass = 0; iPass < 2; iPass++) {
    318           int end = start[2*iPass+1];
    319           for (i = start[2*iPass]; i < end; i++) {
    320                iRow = index[i];
    321                double value = infeas[iRow];
    322                if (value > tolerance) {
    323                     //#define OUT_EQ
    324 #ifdef OUT_EQ
    325                     {
    326                          int iSequence = pivotVariable[iRow];
    327                          if (upper[iSequence] == lower[iSequence])
    328                               value *= 2.0;
    329                     }
    330 #endif
    331                     double weight = CoinMin(weights_[iRow], 1.0e50);
    332                     //largestWeight = CoinMax(largestWeight,weight);
    333                     //smallestWeight = CoinMin(smallestWeight,weight);
    334                     //double dubious = dubiousWeights_[iRow];
    335                     //weight *= dubious;
    336                     //if (value>2.0*largest*weight||(value>0.5*largest*weight&&value*largestWeight>dubious*largest*weight)) {
    337                     if (value > largest * weight) {
    338                          // make last pivot row last resort choice
    339                          if (iRow == lastPivotRow) {
    340                               if (value * 1.0e-10 < largest * weight)
    341                                    continue;
    342                               else
    343                                    value *= 1.0e-10;
    344                          }
    345                          int iSequence = pivotVariable[iRow];
    346                          if (!model_->flagged(iSequence)) {
    347                               //#define CLP_DEBUG 3
    348 #ifdef CLP_DEBUG
    349                               double value2 = 0.0;
    350                               if (solution[iSequence] > upper[iSequence] + tolerance)
    351                                    value2 = solution[iSequence] - upper[iSequence];
    352                               else if (solution[iSequence] < lower[iSequence] - tolerance)
    353                                    value2 = solution[iSequence] - lower[iSequence];
    354                               assert(fabs(value2 * value2 - infeas[iRow]) < 1.0e-8 * CoinMin(value2 * value2, infeas[iRow]));
    355 #endif
    356                               if (solution[iSequence] > upper[iSequence] + tolerance ||
    357                                         solution[iSequence] < lower[iSequence] - tolerance) {
    358                                    chosenRow = iRow;
    359                                    largest = value / weight;
    360                                    //largestWeight = dubious;
    361                               }
    362                          } else {
    363                               // just to make sure we don't exit before got something
    364                               numberWanted++;
    365                          }
    366                     }
    367                     numberWanted--;
    368                     if (!numberWanted)
    369                          break;
    370                }
    371           }
    372           if (!numberWanted)
    373                break;
    374      }
    375      //printf("smallest %g largest %g\n",smallestWeight,largestWeight);
    376      if (chosenRow < 0 && toleranceChanged) {
    377           // won't line up with checkPrimalSolution - do again
    378           double saveError = model_->largestDualError();
    379           model_->setLargestDualError(0.0);
    380           // can't loop
    381           chosenRow = pivotRow();
    382           model_->setLargestDualError(saveError);
    383      }
    384      if (chosenRow < 0 && lastPivotRow < 0) {
    385           int nLeft = 0;
    386           for (int i = 0; i < number; i++) {
    387                int iRow = index[i];
    388                if (fabs(infeas[iRow]) > 1.0e-50) {
    389                     index[nLeft++] = iRow;
    390                } else {
    391                     infeas[iRow] = 0.0;
    392                }
    393           }
    394           infeasible_->setNumElements(nLeft);
    395           model_->setNumberPrimalInfeasibilities(nLeft);
    396      }
    397      return chosenRow;
     362        }
     363        numberWanted--;
     364        if (!numberWanted)
     365          break;
     366      }
     367    }
     368    if (!numberWanted)
     369      break;
     370  }
     371  //printf("smallest %g largest %g\n",smallestWeight,largestWeight);
     372  if (chosenRow < 0 && toleranceChanged) {
     373    // won't line up with checkPrimalSolution - do again
     374    double saveError = model_->largestDualError();
     375    model_->setLargestDualError(0.0);
     376    // can't loop
     377    chosenRow = pivotRow();
     378    model_->setLargestDualError(saveError);
     379  }
     380  if (chosenRow < 0 && lastPivotRow < 0) {
     381    int nLeft = 0;
     382    for (int i = 0; i < number; i++) {
     383      int iRow = index[i];
     384      if (fabs(infeas[iRow]) > 1.0e-50) {
     385        index[nLeft++] = iRow;
     386      } else {
     387        infeas[iRow] = 0.0;
     388      }
     389    }
     390    infeasible_->setNumElements(nLeft);
     391    model_->setNumberPrimalInfeasibilities(nLeft);
     392  }
     393  return chosenRow;
    398394}
    399395#if 0
     
    407403   Also does FT update */
    408404double
    409 ClpDualRowSteepest::updateWeights(CoinIndexedVector * input,
    410                                   CoinIndexedVector * spare,
    411                                   CoinIndexedVector * spare2,
    412                                   CoinIndexedVector * updatedColumn)
    413 {
    414      //#define CLP_DEBUG 3
    415 #if CLP_DEBUG>2
    416      // Very expensive debug
    417      {
    418           int numberRows = model_->numberRows();
    419           CoinIndexedVector * temp = new CoinIndexedVector();
    420           temp->reserve(numberRows +
    421                         model_->factorization()->maximumPivots());
    422           double * array = alternateWeights_->denseVector();
    423           int * which = alternateWeights_->getIndices();
    424           alternateWeights_->clear();
    425           int i;
    426           for (i = 0; i < numberRows; i++) {
    427                double value = 0.0;
    428                array[i] = 1.0;
    429                which[0] = i;
    430                alternateWeights_->setNumElements(1);
    431                model_->factorization()->updateColumnTranspose(temp,
    432                          alternateWeights_);
    433                int number = alternateWeights_->getNumElements();
    434                int j;
    435                for (j = 0; j < number; j++) {
    436                     int iRow = which[j];
    437                     value += array[iRow] * array[iRow];
    438                     array[iRow] = 0.0;
    439                }
    440                alternateWeights_->setNumElements(0);
    441                double w = CoinMax(weights_[i], value) * .1;
    442                if (fabs(weights_[i] - value) > w) {
    443                     printf("%d old %g, true %g\n", i, weights_[i], value);
    444                     weights_[i] = value; // to reduce printout
    445                }
    446                //else
    447                //printf("%d matches %g\n",i,value);
    448           }
    449           delete temp;
    450      }
    451 #endif
    452      assert (input->packedMode());
    453      if (!updatedColumn->packedMode()) {
    454           // I think this means empty
     405ClpDualRowSteepest::updateWeights(CoinIndexedVector *input,
     406  CoinIndexedVector *spare,
     407  CoinIndexedVector *spare2,
     408  CoinIndexedVector *updatedColumn)
     409{
     410  //#define CLP_DEBUG 3
     411#if CLP_DEBUG > 2
     412  // Very expensive debug
     413  {
     414    int numberRows = model_->numberRows();
     415    CoinIndexedVector *temp = new CoinIndexedVector();
     416    temp->reserve(numberRows + model_->factorization()->maximumPivots());
     417    double *array = alternateWeights_->denseVector();
     418    int *which = alternateWeights_->getIndices();
     419    alternateWeights_->clear();
     420    int i;
     421    for (i = 0; i < numberRows; i++) {
     422      double value = 0.0;
     423      array[i] = 1.0;
     424      which[0] = i;
     425      alternateWeights_->setNumElements(1);
     426      model_->factorization()->updateColumnTranspose(temp,
     427        alternateWeights_);
     428      int number = alternateWeights_->getNumElements();
     429      int j;
     430      for (j = 0; j < number; j++) {
     431        int iRow = which[j];
     432        value += array[iRow] * array[iRow];
     433        array[iRow] = 0.0;
     434      }
     435      alternateWeights_->setNumElements(0);
     436      double w = CoinMax(weights_[i], value) * .1;
     437      if (fabs(weights_[i] - value) > w) {
     438        printf("%d old %g, true %g\n", i, weights_[i], value);
     439        weights_[i] = value; // to reduce printout
     440      }
     441      //else
     442      //printf("%d matches %g\n",i,value);
     443    }
     444    delete temp;
     445  }
     446#endif
     447  assert(input->packedMode());
     448  if (!updatedColumn->packedMode()) {
     449    // I think this means empty
    455450#ifdef COIN_DEVELOP
    456           printf("updatedColumn not packed mode ClpDualRowSteepest::updateWeights\n");
    457 #endif
    458           return 0.0;
    459      }
    460      double alpha = 0.0;
    461      if (!model_->factorization()->networkBasis()) {
    462           // clear other region
    463           alternateWeights_->clear();
    464           double norm = 0.0;
    465           int i;
    466           double * work = input->denseVector();
    467           int numberNonZero = input->getNumElements();
    468           int * which = input->getIndices();
    469           double * work2 = spare->denseVector();
    470           int * which2 = spare->getIndices();
    471           // ftran
    472           //permute and move indices into index array
    473           //also compute norm
    474           //int *regionIndex = alternateWeights_->getIndices (  );
    475           const int *permute = model_->factorization()->permute();
    476           //double * region = alternateWeights_->denseVector();
    477           if (permute) {
    478                for ( i = 0; i < numberNonZero; i ++ ) {
    479                     int iRow = which[i];
    480                     double value = work[i];
    481                     norm += value * value;
    482                     iRow = permute[iRow];
    483                     work2[iRow] = value;
    484                     which2[i] = iRow;
    485                }
    486           } else {
    487                for ( i = 0; i < numberNonZero; i ++ ) {
    488                     int iRow = which[i];
    489                     double value = work[i];
    490                     norm += value * value;
    491                     //iRow = permute[iRow];
    492                     work2[iRow] = value;
    493                     which2[i] = iRow;
    494                }
    495           }
    496           spare->setNumElements ( numberNonZero );
    497           // Do FT update
     451    printf("updatedColumn not packed mode ClpDualRowSteepest::updateWeights\n");
     452#endif
     453    return 0.0;
     454  }
     455  double alpha = 0.0;
     456  if (!model_->factorization()->networkBasis()) {
     457    // clear other region
     458    alternateWeights_->clear();
     459    double norm = 0.0;
     460    int i;
     461    double *work = input->denseVector();
     462    int numberNonZero = input->getNumElements();
     463    int *which = input->getIndices();
     464    double *work2 = spare->denseVector();
     465    int *which2 = spare->getIndices();
     466    // ftran
     467    //permute and move indices into index array
     468    //also compute norm
     469    //int *regionIndex = alternateWeights_->getIndices (  );
     470    const int *permute = model_->factorization()->permute();
     471    //double * region = alternateWeights_->denseVector();
     472    if (permute) {
     473      for (i = 0; i < numberNonZero; i++) {
     474        int iRow = which[i];
     475        double value = work[i];
     476        norm += value * value;
     477        iRow = permute[iRow];
     478        work2[iRow] = value;
     479        which2[i] = iRow;
     480      }
     481    } else {
     482      for (i = 0; i < numberNonZero; i++) {
     483        int iRow = which[i];
     484        double value = work[i];
     485        norm += value * value;
     486        //iRow = permute[iRow];
     487        work2[iRow] = value;
     488        which2[i] = iRow;
     489      }
     490    }
     491    spare->setNumElements(numberNonZero);
     492    // Do FT update
    498493#if 0
    499494          ft_count_in += updatedColumn->getNumElements();
    500495          up_count_in += spare->getNumElements();
    501496#endif
    502           if (permute || true) {
    503 #if CLP_DEBUG>2
    504                printf("REGION before %d els\n", spare->getNumElements());
    505                spare->print();
    506 #endif
    507                model_->factorization()->updateTwoColumnsFT(spare2, updatedColumn,
    508                          spare, permute != NULL);
    509 #if CLP_DEBUG>2
    510                printf("REGION after %d els\n", spare->getNumElements());
    511                spare->print();
    512 #endif
    513           } else {
    514                // Leave as old way
    515                model_->factorization()->updateColumnFT(spare2, updatedColumn);
    516                model_->factorization()->updateColumn(spare2, spare, false);
    517           }
     497    if (permute || true) {
     498#if CLP_DEBUG > 2
     499      printf("REGION before %d els\n", spare->getNumElements());
     500      spare->print();
     501#endif
     502      model_->factorization()->updateTwoColumnsFT(spare2, updatedColumn,
     503        spare, permute != NULL);
     504#if CLP_DEBUG > 2
     505      printf("REGION after %d els\n", spare->getNumElements());
     506      spare->print();
     507#endif
     508    } else {
     509      // Leave as old way
     510      model_->factorization()->updateColumnFT(spare2, updatedColumn);
     511      model_->factorization()->updateColumn(spare2, spare, false);
     512    }
    518513#undef CLP_DEBUG
    519514#if 0
     
    525520                      ft_count_in, up_count_in);
    526521#endif
    527           numberNonZero = spare->getNumElements();
    528           // alternateWeights_ should still be empty
    529           int pivotRow = model_->pivotRow();
     522    numberNonZero = spare->getNumElements();
     523    // alternateWeights_ should still be empty
     524    int pivotRow = model_->pivotRow();
    530525#ifdef CLP_DEBUG
    531           if ( model_->logLevel (  ) > 4  &&
    532                     fabs(norm - weights_[pivotRow]) > 1.0e-3 * (1.0 + norm))
    533                printf("on row %d, true weight %g, old %g\n",
    534                       pivotRow, sqrt(norm), sqrt(weights_[pivotRow]));
    535 #endif
    536           // could re-initialize here (could be expensive)
    537           norm /= model_->alpha() * model_->alpha();
    538           assert(model_->alpha());
    539           assert(norm);
    540           // pivot element
    541           alpha = 0.0;
    542           double multiplier = 2.0 / model_->alpha();
    543           // look at updated column
    544           work = updatedColumn->denseVector();
    545           numberNonZero = updatedColumn->getNumElements();
    546           which = updatedColumn->getIndices();
     526    if (model_->logLevel() > 4 && fabs(norm - weights_[pivotRow]) > 1.0e-3 * (1.0 + norm))
     527      printf("on row %d, true weight %g, old %g\n",
     528        pivotRow, sqrt(norm), sqrt(weights_[pivotRow]));
     529#endif
     530    // could re-initialize here (could be expensive)
     531    norm /= model_->alpha() * model_->alpha();
     532    assert(model_->alpha());
     533    assert(norm);
     534    // pivot element
     535    alpha = 0.0;
     536    double multiplier = 2.0 / model_->alpha();
     537    // look at updated column
     538    work = updatedColumn->denseVector();
     539    numberNonZero = updatedColumn->getNumElements();
     540    which = updatedColumn->getIndices();
    547541
    548           int nSave = 0;
    549           double * work3 = alternateWeights_->denseVector();
    550           int * which3 = alternateWeights_->getIndices();
    551           const int * pivotColumn = model_->factorization()->pivotColumn();
    552           for (i = 0; i < numberNonZero; i++) {
    553                int iRow = which[i];
    554                double theta = work[i];
    555                if (iRow == pivotRow)
    556                     alpha = theta;
    557                double devex = weights_[iRow];
    558                work3[nSave] = devex; // save old
    559                which3[nSave++] = iRow;
    560                // transform to match spare
    561                int jRow = permute ? pivotColumn[iRow] : iRow;
    562                double value = work2[jRow];
    563                devex += theta * (theta * norm + value * multiplier);
    564                if (devex < DEVEX_TRY_NORM)
    565                     devex = DEVEX_TRY_NORM;
    566                weights_[iRow] = devex;
    567           }
    568           alternateWeights_->setPackedMode(true);
    569           alternateWeights_->setNumElements(nSave);
    570           if (norm < DEVEX_TRY_NORM)
    571                norm = DEVEX_TRY_NORM;
    572           // Try this to make less likely will happen again and stop cycling
    573           //norm *= 1.02;
    574           weights_[pivotRow] = norm;
    575           spare->clear();
     542    int nSave = 0;
     543    double *work3 = alternateWeights_->denseVector();
     544    int *which3 = alternateWeights_->getIndices();
     545    const int *pivotColumn = model_->factorization()->pivotColumn();
     546    for (i = 0; i < numberNonZero; i++) {
     547      int iRow = which[i];
     548      double theta = work[i];
     549      if (iRow == pivotRow)
     550        alpha = theta;
     551      double devex = weights_[iRow];
     552      work3[nSave] = devex; // save old
     553      which3[nSave++] = iRow;
     554      // transform to match spare
     555      int jRow = permute ? pivotColumn[iRow] : iRow;
     556      double value = work2[jRow];
     557      devex += theta * (theta * norm + value * multiplier);
     558      if (devex < DEVEX_TRY_NORM)
     559        devex = DEVEX_TRY_NORM;
     560      weights_[iRow] = devex;
     561    }
     562    alternateWeights_->setPackedMode(true);
     563    alternateWeights_->setNumElements(nSave);
     564    if (norm < DEVEX_TRY_NORM)
     565      norm = DEVEX_TRY_NORM;
     566    // Try this to make less likely will happen again and stop cycling
     567    //norm *= 1.02;
     568    weights_[pivotRow] = norm;
     569    spare->clear();
    576570#ifdef CLP_DEBUG
    577           spare->checkClear();
    578 #endif
    579      } else {
    580           // Do FT update
    581           model_->factorization()->updateColumnFT(spare, updatedColumn);
    582           // clear other region
    583           alternateWeights_->clear();
    584           double norm = 0.0;
    585           int i;
    586           double * work = input->denseVector();
    587           int number = input->getNumElements();
    588           int * which = input->getIndices();
    589           double * work2 = spare->denseVector();
    590           int * which2 = spare->getIndices();
    591           for (i = 0; i < number; i++) {
    592                int iRow = which[i];
    593                double value = work[i];
    594                norm += value * value;
    595                work2[iRow] = value;
    596                which2[i] = iRow;
    597           }
    598           spare->setNumElements(number);
    599           // ftran
     571    spare->checkClear();
     572#endif
     573  } else {
     574    // Do FT update
     575    model_->factorization()->updateColumnFT(spare, updatedColumn);
     576    // clear other region
     577    alternateWeights_->clear();
     578    double norm = 0.0;
     579    int i;
     580    double *work = input->denseVector();
     581    int number = input->getNumElements();
     582    int *which = input->getIndices();
     583    double *work2 = spare->denseVector();
     584    int *which2 = spare->getIndices();
     585    for (i = 0; i < number; i++) {
     586      int iRow = which[i];
     587      double value = work[i];
     588      norm += value * value;
     589      work2[iRow] = value;
     590      which2[i] = iRow;
     591    }
     592    spare->setNumElements(number);
     593    // ftran
    600594#ifndef NDEBUG
    601           alternateWeights_->checkClear();
    602 #endif
    603           model_->factorization()->updateColumn(alternateWeights_, spare);
    604           // alternateWeights_ should still be empty
     595    alternateWeights_->checkClear();
     596#endif
     597    model_->factorization()->updateColumn(alternateWeights_, spare);
     598    // alternateWeights_ should still be empty
    605599#ifndef NDEBUG
    606           alternateWeights_->checkClear();
    607 #endif
    608           int pivotRow = model_->pivotRow();
     600    alternateWeights_->checkClear();
     601#endif
     602    int pivotRow = model_->pivotRow();
    609603#ifdef CLP_DEBUG
    610           if ( model_->logLevel (  ) > 4  &&
    611                     fabs(norm - weights_[pivotRow]) > 1.0e-3 * (1.0 + norm))
    612                printf("on row %d, true weight %g, old %g\n",
    613                       pivotRow, sqrt(norm), sqrt(weights_[pivotRow]));
    614 #endif
    615           // could re-initialize here (could be expensive)
    616           norm /= model_->alpha() * model_->alpha();
     604    if (model_->logLevel() > 4 && fabs(norm - weights_[pivotRow]) > 1.0e-3 * (1.0 + norm))
     605      printf("on row %d, true weight %g, old %g\n",
     606        pivotRow, sqrt(norm), sqrt(weights_[pivotRow]));
     607#endif
     608    // could re-initialize here (could be expensive)
     609    norm /= model_->alpha() * model_->alpha();
    617610
    618           assert(norm);
    619           //if (norm < DEVEX_TRY_NORM)
    620           //norm = DEVEX_TRY_NORM;
    621           // pivot element
    622           alpha = 0.0;
    623           double multiplier = 2.0 / model_->alpha();
    624           // look at updated column
    625           work = updatedColumn->denseVector();
    626           number = updatedColumn->getNumElements();
    627           which = updatedColumn->getIndices();
     611    assert(norm);
     612    //if (norm < DEVEX_TRY_NORM)
     613    //norm = DEVEX_TRY_NORM;
     614    // pivot element
     615    alpha = 0.0;
     616    double multiplier = 2.0 / model_->alpha();
     617    // look at updated column
     618    work = updatedColumn->denseVector();
     619    number = updatedColumn->getNumElements();
     620    which = updatedColumn->getIndices();
    628621
    629           int nSave = 0;
    630           double * work3 = alternateWeights_->denseVector();
    631           int * which3 = alternateWeights_->getIndices();
    632           for (i = 0; i < number; i++) {
    633                int iRow = which[i];
    634                double theta = work[i];
    635                if (iRow == pivotRow)
    636                     alpha = theta;
    637                double devex = weights_[iRow];
    638                work3[nSave] = devex; // save old
    639                which3[nSave++] = iRow;
    640                double value = work2[iRow];
    641                devex += theta * (theta * norm + value * multiplier);
    642                if (devex < DEVEX_TRY_NORM)
    643                     devex = DEVEX_TRY_NORM;
    644                weights_[iRow] = devex;
    645           }
    646           if (!alpha) {
    647                // error - but carry on
    648                alpha = 1.0e-50;
    649           }
    650           alternateWeights_->setPackedMode(true);
    651           alternateWeights_->setNumElements(nSave);
    652           if (norm < DEVEX_TRY_NORM)
    653                norm = DEVEX_TRY_NORM;
    654           weights_[pivotRow] = norm;
    655           spare->clear();
    656      }
     622    int nSave = 0;
     623    double *work3 = alternateWeights_->denseVector();
     624    int *which3 = alternateWeights_->getIndices();
     625    for (i = 0; i < number; i++) {
     626      int iRow = which[i];
     627      double theta = work[i];
     628      if (iRow == pivotRow)
     629        alpha = theta;
     630      double devex = weights_[iRow];
     631      work3[nSave] = devex; // save old
     632      which3[nSave++] = iRow;
     633      double value = work2[iRow];
     634      devex += theta * (theta * norm + value * multiplier);
     635      if (devex < DEVEX_TRY_NORM)
     636        devex = DEVEX_TRY_NORM;
     637      weights_[iRow] = devex;
     638    }
     639    if (!alpha) {
     640      // error - but carry on
     641      alpha = 1.0e-50;
     642    }
     643    alternateWeights_->setPackedMode(true);
     644    alternateWeights_->setNumElements(nSave);
     645    if (norm < DEVEX_TRY_NORM)
     646      norm = DEVEX_TRY_NORM;
     647    weights_[pivotRow] = norm;
     648    spare->clear();
     649  }
    657650#ifdef CLP_DEBUG
    658      spare->checkClear();
    659 #endif
    660      return alpha;
     651  spare->checkClear();
     652#endif
     653  return alpha;
    661654}
    662655
     
    665658   Computes change in objective function
    666659*/
    667 void
    668 ClpDualRowSteepest::updatePrimalSolution(
    669      CoinIndexedVector * primalUpdate,
    670      double primalRatio,
    671      double & objectiveChange)
    672 {
    673      double * COIN_RESTRICT work = primalUpdate->denseVector();
    674      int number = primalUpdate->getNumElements();
    675      int * COIN_RESTRICT which = primalUpdate->getIndices();
    676      int i;
    677      double changeObj = 0.0;
    678      double tolerance = model_->currentPrimalTolerance();
    679      const int * COIN_RESTRICT pivotVariable = model_->pivotVariable();
    680      double * COIN_RESTRICT infeas = infeasible_->denseVector();
    681      double * COIN_RESTRICT solution = model_->solutionRegion();
    682      const double * COIN_RESTRICT costModel = model_->costRegion();
    683      const double * COIN_RESTRICT lowerModel = model_->lowerRegion();
    684      const double * COIN_RESTRICT upperModel = model_->upperRegion();
     660void ClpDualRowSteepest::updatePrimalSolution(
     661  CoinIndexedVector *primalUpdate,
     662  double primalRatio,
     663  double &objectiveChange)
     664{
     665  double *COIN_RESTRICT work = primalUpdate->denseVector();
     666  int number = primalUpdate->getNumElements();
     667  int *COIN_RESTRICT which = primalUpdate->getIndices();
     668  int i;
     669  double changeObj = 0.0;
     670  double tolerance = model_->currentPrimalTolerance();
     671  const int *COIN_RESTRICT pivotVariable = model_->pivotVariable();
     672  double *COIN_RESTRICT infeas = infeasible_->denseVector();
     673  double *COIN_RESTRICT solution = model_->solutionRegion();
     674  const double *COIN_RESTRICT costModel = model_->costRegion();
     675  const double *COIN_RESTRICT lowerModel = model_->lowerRegion();
     676  const double *COIN_RESTRICT upperModel = model_->upperRegion();
    685677#ifdef CLP_DUAL_COLUMN_MULTIPLIER
    686      int numberColumns = model_->numberColumns();
    687 #endif
    688      if (primalUpdate->packedMode()) {
    689           for (i = 0; i < number; i++) {
    690                int iRow = which[i];
    691                int iPivot = pivotVariable[iRow];
    692                double value = solution[iPivot];
    693                double cost = costModel[iPivot];
    694                double change = primalRatio * work[i];
    695                work[i] = 0.0;
    696                value -= change;
    697                changeObj -= change * cost;
    698                double lower = lowerModel[iPivot];
    699                double upper = upperModel[iPivot];
    700                solution[iPivot] = value;
    701                if (value < lower - tolerance) {
    702                     value -= lower;
    703                     value *= value;
     678  int numberColumns = model_->numberColumns();
     679#endif
     680  if (primalUpdate->packedMode()) {
     681    for (i = 0; i < number; i++) {
     682      int iRow = which[i];
     683      int iPivot = pivotVariable[iRow];
     684      double value = solution[iPivot];
     685      double cost = costModel[iPivot];
     686      double change = primalRatio * work[i];
     687      work[i] = 0.0;
     688      value -= change;
     689      changeObj -= change * cost;
     690      double lower = lowerModel[iPivot];
     691      double upper = upperModel[iPivot];
     692      solution[iPivot] = value;
     693      if (value < lower - tolerance) {
     694        value -= lower;
     695        value *= value;
    704696#ifdef CLP_DUAL_COLUMN_MULTIPLIER
    705                     if (iPivot < numberColumns)
    706                          value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
     697        if (iPivot < numberColumns)
     698          value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
    707699#endif
    708700#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
    709                     if (lower == upper)
    710                          value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
    711 #endif
    712                     // store square in list
    713                     if (infeas[iRow])
    714                          infeas[iRow] = value; // already there
    715                     else
    716                          infeasible_->quickAdd(iRow, value);
    717                } else if (value > upper + tolerance) {
    718                     value -= upper;
    719                     value *= value;
     701        if (lower == upper)
     702          value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
     703#endif
     704        // store square in list
     705        if (infeas[iRow])
     706          infeas[iRow] = value; // already there
     707        else
     708          infeasible_->quickAdd(iRow, value);
     709      } else if (value > upper + tolerance) {
     710        value -= upper;
     711        value *= value;
    720712#ifdef CLP_DUAL_COLUMN_MULTIPLIER
    721                     if (iPivot < numberColumns)
    722                          value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
     713        if (iPivot < numberColumns)
     714          value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
    723715#endif
    724716#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
    725                     if (lower == upper)
    726                          value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
    727 #endif
    728                     // store square in list
    729                     if (infeas[iRow])
    730                          infeas[iRow] = value; // already there
    731                     else
    732                          infeasible_->quickAdd(iRow, value);
    733                } else {
    734                     // feasible - was it infeasible - if so set tiny
    735                     if (infeas[iRow])
    736                          infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
    737                }
    738           }
    739      } else {
    740           for (i = 0; i < number; i++) {
    741                int iRow = which[i];
    742                int iPivot = pivotVariable[iRow];
    743                double value = solution[iPivot];
    744                double cost = costModel[iPivot];
    745                double change = primalRatio * work[iRow];
    746                value -= change;
    747                changeObj -= change * cost;
    748                double lower = lowerModel[iPivot];
    749                double upper = upperModel[iPivot];
    750                solution[iPivot] = value;
    751                if (value < lower - tolerance) {
    752                     value -= lower;
    753                     value *= value;
     717        if (lower == upper)
     718          value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
     719#endif
     720        // store square in list
     721        if (infeas[iRow])
     722          infeas[iRow] = value; // already there
     723        else
     724          infeasible_->quickAdd(iRow, value);
     725      } else {
     726        // feasible - was it infeasible - if so set tiny
     727        if (infeas[iRow])
     728          infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
     729      }
     730    }
     731  } else {
     732    for (i = 0; i < number; i++) {
     733      int iRow = which[i];
     734      int iPivot = pivotVariable[iRow];
     735      double value = solution[iPivot];
     736      double cost = costModel[iPivot];
     737      double change = primalRatio * work[iRow];
     738      value -= change;
     739      changeObj -= change * cost;
     740      double lower = lowerModel[iPivot];
     741      double upper = upperModel[iPivot];
     742      solution[iPivot] = value;
     743      if (value < lower - tolerance) {
     744        value -= lower;
     745        value *= value;
    754746#ifdef CLP_DUAL_COLUMN_MULTIPLIER
    755                     if (iPivot < numberColumns)
    756                          value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
     747        if (iPivot < numberColumns)
     748          value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
    757749#endif
    758750#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
    759                     if (lower == upper)
    760                          value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
    761 #endif
    762                     // store square in list
    763                     if (infeas[iRow])
    764                          infeas[iRow] = value; // already there
    765                     else
    766                          infeasible_->quickAdd(iRow, value);
    767                } else if (value > upper + tolerance) {
    768                     value -= upper;
    769                     value *= value;
     751        if (lower == upper)
     752          value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
     753#endif
     754        // store square in list
     755        if (infeas[iRow])
     756          infeas[iRow] = value; // already there
     757        else
     758          infeasible_->quickAdd(iRow, value);
     759      } else if (value > upper + tolerance) {
     760        value -= upper;
     761        value *= value;
    770762#ifdef CLP_DUAL_COLUMN_MULTIPLIER
    771                     if (iPivot < numberColumns)
    772                          value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
     763        if (iPivot < numberColumns)
     764          value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
    773765#endif
    774766#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
    775                     if (lower == upper)
    776                          value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
    777 #endif
    778                     // store square in list
    779                     if (infeas[iRow])
    780                          infeas[iRow] = value; // already there
    781                     else
    782                          infeasible_->quickAdd(iRow, value);
    783                } else {
    784                     // feasible - was it infeasible - if so set tiny
    785                     if (infeas[iRow])
    786                          infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
    787                }
    788                work[iRow] = 0.0;
    789           }
    790      }
    791      // Do pivot row
    792      {
    793           int iRow = model_->pivotRow();
    794           // feasible - was it infeasible - if so set tiny
    795           //assert (infeas[iRow]);
    796           if (infeas[iRow])
    797                infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
    798      }
    799      primalUpdate->setNumElements(0);
    800      objectiveChange += changeObj;
     767        if (lower == upper)
     768          value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
     769#endif
     770        // store square in list
     771        if (infeas[iRow])
     772          infeas[iRow] = value; // already there
     773        else
     774          infeasible_->quickAdd(iRow, value);
     775      } else {
     776        // feasible - was it infeasible - if so set tiny
     777        if (infeas[iRow])
     778          infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
     779      }
     780      work[iRow] = 0.0;
     781    }
     782  }
     783  // Do pivot row
     784  {
     785    int iRow = model_->pivotRow();
     786    // feasible - was it infeasible - if so set tiny
     787    //assert (infeas[iRow]);
     788    if (infeas[iRow])
     789      infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
     790  }
     791  primalUpdate->setNumElements(0);
     792  objectiveChange += changeObj;
    801793}
    802794/* Saves any weights round factorization as pivot rows may change
     
    809801   7) for strong branching - initialize full weights , infeasibilities
    810802*/
    811 void
    812 ClpDualRowSteepest::saveWeights(ClpSimplex * model, int mode)
    813 {
    814      // alternateWeights_ is defined as indexed but is treated oddly
    815      model_ = model;
    816      int numberRows = model_->numberRows();
    817      int numberColumns = model_->numberColumns();
    818      const int * pivotVariable = model_->pivotVariable();
    819      int i;
    820      if (mode == 1) {
    821           if(weights_) {
    822                // Check if size has changed
    823                if (infeasible_->capacity() == numberRows) {
    824                     alternateWeights_->clear();
    825                     // change from row numbers to sequence numbers
    826                     int * which = alternateWeights_->getIndices();
    827                     for (i = 0; i < numberRows; i++) {
    828                          int iPivot = pivotVariable[i];
    829                          which[i] = iPivot;
    830                     }
    831                     state_ = 1;
    832                     // clear marker on savedWeights_
    833                     assert (savedWeights_);
    834                     if (savedWeights_->packedMode()) {
    835                       savedWeights_->setPackedMode(false);
    836                       savedWeights_->setNumElements(0);
    837                     }
    838                } else {
    839                     // size has changed - clear everything
    840                     delete [] weights_;
    841                     weights_ = NULL;
    842                     delete [] dubiousWeights_;
    843                     dubiousWeights_ = NULL;
    844                     delete infeasible_;
    845                     infeasible_ = NULL;
    846                     delete alternateWeights_;
    847                     alternateWeights_ = NULL;
    848                     delete savedWeights_;
    849                     savedWeights_ = NULL;
    850                     state_ = -1;
    851                }
     803void ClpDualRowSteepest::saveWeights(ClpSimplex *model, int mode)
     804{
     805  // alternateWeights_ is defined as indexed but is treated oddly
     806  model_ = model;
     807  int numberRows = model_->numberRows();
     808  int numberColumns = model_->numberColumns();
     809  const int *pivotVariable = model_->pivotVariable();
     810  int i;
     811  if (mode == 1) {
     812    if (weights_) {
     813      // Check if size has changed
     814      if (infeasible_->capacity() == numberRows) {
     815        alternateWeights_->clear();
     816        // change from row numbers to sequence numbers
     817        int *which = alternateWeights_->getIndices();
     818        for (i = 0; i < numberRows; i++) {
     819          int iPivot = pivotVariable[i];
     820          which[i] = iPivot;
     821        }
     822        state_ = 1;
     823        // clear marker on savedWeights_
     824        assert(savedWeights_);
     825        if (savedWeights_->packedMode()) {
     826          savedWeights_->setPackedMode(false);
     827          savedWeights_->setNumElements(0);
     828        }
     829      } else {
     830        // size has changed - clear everything
     831        delete[] weights_;
     832        weights_ = NULL;
     833        delete[] dubiousWeights_;
     834        dubiousWeights_ = NULL;
     835        delete infeasible_;
     836        infeasible_ = NULL;
     837        delete alternateWeights_;
     838        alternateWeights_ = NULL;
     839        delete savedWeights_;
     840        savedWeights_ = NULL;
     841        state_ = -1;
     842      }
     843    }
     844  } else if (mode == 2 || mode == 4 || mode >= 5) {
     845    // restore
     846    if (!weights_ || state_ == -1 || mode == 5 || mode == 7) {
     847      // initialize weights
     848      delete[] weights_;
     849      delete alternateWeights_;
     850      weights_ = new double[numberRows];
     851      // initialize to 1.0 (can we do better?)
     852      for (i = 0; i < numberRows; i++) {
     853        weights_[i] = 1.0;
     854      }
     855      alternateWeights_ = new CoinIndexedVector();
     856      // enough space so can use it for factorization
     857      alternateWeights_->reserve(numberRows + model_->factorization()->maximumPivots());
     858      if (mode_ != 1 || mode == 5) {
     859      } else {
     860        CoinIndexedVector *temp = new CoinIndexedVector();
     861        temp->reserve(numberRows + model_->factorization()->maximumPivots());
     862        double *array = alternateWeights_->denseVector();
     863        int *which = alternateWeights_->getIndices();
     864        int firstRow = 0;
     865        int lastRow = numberRows;
     866        if (mode == 7) {
     867          // use info passed in
     868          firstRow = model->spareIntArray_[0];
     869          lastRow = model->spareIntArray_[1];
     870        }
     871        for (i = firstRow; i < lastRow; i++) {
     872          double value = 0.0;
     873          array[0] = 1.0;
     874          which[0] = i;
     875          alternateWeights_->setNumElements(1);
     876          alternateWeights_->setPackedMode(true);
     877          model_->factorization()->updateColumnTranspose(temp,
     878            alternateWeights_);
     879          int number = alternateWeights_->getNumElements();
     880          int j;
     881          for (j = 0; j < number; j++) {
     882            value += array[j] * array[j];
     883            array[j] = 0.0;
    852884          }
    853      } else if (mode == 2 || mode == 4 || mode >= 5) {
    854           // restore
    855           if (!weights_ || state_ == -1 || mode == 5 || mode == 7) {
    856                // initialize weights
    857                delete [] weights_;
    858                delete alternateWeights_;
    859                weights_ = new double[numberRows];
    860                // initialize to 1.0 (can we do better?)
    861                for (i = 0; i < numberRows; i++) {
    862                  weights_[i] = 1.0;
    863                }
    864                alternateWeights_ = new CoinIndexedVector();
    865                // enough space so can use it for factorization
    866                alternateWeights_->reserve(numberRows +
    867                                           model_->factorization()->maximumPivots());
    868                if (mode_ != 1 || mode == 5) {
    869                } else {
    870                     CoinIndexedVector * temp = new CoinIndexedVector();
    871                     temp->reserve(numberRows +
    872                                   model_->factorization()->maximumPivots());
    873                     double * array = alternateWeights_->denseVector();
    874                     int * which = alternateWeights_->getIndices();
    875                     int firstRow=0;
    876                     int lastRow=numberRows;
    877                     if (mode==7) {
    878                       // use info passed in
    879                       firstRow=model->spareIntArray_[0];
    880                       lastRow=model->spareIntArray_[1];
    881                     }
    882                     for (i = firstRow; i < lastRow; i++) {
    883                          double value = 0.0;
    884                          array[0] = 1.0;
    885                          which[0] = i;
    886                          alternateWeights_->setNumElements(1);
    887                          alternateWeights_->setPackedMode(true);
    888                          model_->factorization()->updateColumnTranspose(temp,
    889                                    alternateWeights_);
    890                          int number = alternateWeights_->getNumElements();
    891                          int j;
    892                          for (j = 0; j < number; j++) {
    893                               value += array[j] * array[j];
    894                               array[j] = 0.0;
    895                          }
    896                          alternateWeights_->setNumElements(0);
    897                          weights_[i] = value;
    898                     }
    899                     delete temp;
    900                }
    901                // create saved weights (not really indexedvector)
    902                savedWeights_ = new CoinIndexedVector();
    903                savedWeights_->reserve(numberRows);
    904                for (int i=0;i<model_->numberRows();i++)
    905                  savedWeights_->denseVector()[i]=1.0;
     885          alternateWeights_->setNumElements(0);
     886          weights_[i] = value;
     887        }
     888        delete temp;
     889      }
     890      // create saved weights (not really indexedvector)
     891      savedWeights_ = new CoinIndexedVector();
     892      savedWeights_->reserve(numberRows);
     893      for (int i = 0; i < model_->numberRows(); i++)
     894        savedWeights_->denseVector()[i] = 1.0;
    906895
    907                double * array = savedWeights_->denseVector();
    908                int * which = savedWeights_->getIndices();
    909                for (int i = 0; i < numberRows; i++) {
    910                     array[i] = weights_[i];
    911                     which[i] = pivotVariable[i];
    912                }
    913                if (mode==7) {
    914                 savedWeights_->setNumElements(numberRows); //flag as special
    915                 savedWeights_->setPackedMode(true);
    916                }
    917           } else if (mode != 6) {
    918                int * which = alternateWeights_->getIndices();
    919                CoinIndexedVector * rowArray3 = model_->rowArray(3);
    920                rowArray3->clear();
    921                int * back = rowArray3->getIndices();
    922                // In case something went wrong
    923                for (i = 0; i < numberRows + numberColumns; i++)
    924                     back[i] = -1;
    925                if (mode != 4) {
    926                     // save
    927                     CoinMemcpyN(which,  numberRows, savedWeights_->getIndices());
    928                     CoinMemcpyN(weights_,       numberRows, savedWeights_->denseVector());
    929                } else {
    930                     // restore
    931                     //memcpy(which,savedWeights_->getIndices(),
    932                     //     numberRows*sizeof(int));
    933                     //memcpy(weights_,savedWeights_->denseVector(),
    934                     //     numberRows*sizeof(double));
    935                     which = savedWeights_->getIndices();
    936                }
    937                // restore (a bit slow - but only every re-factorization)
    938                double * array = savedWeights_->denseVector();
    939                for (i = 0; i < numberRows; i++) {
    940                     int iSeq = which[i];
    941                     back[iSeq] = i;
    942                }
    943                int firstRow=0;
    944                int lastRow=numberRows;
    945                if (mode==7) {
    946                 // use info passed in
    947                  firstRow=model->spareIntArray_[0];
    948                  lastRow=model->spareIntArray_[1];
    949                }
    950                for (i = firstRow; i < lastRow; i++) {
    951                     int iPivot = pivotVariable[i];
    952                     iPivot = back[iPivot];
    953                     if (iPivot >= 0) {
    954                          weights_[i] = array[iPivot];
    955                          if (weights_[i] < DEVEX_TRY_NORM)
    956                               weights_[i] = DEVEX_TRY_NORM; // may need to check more
    957                     } else {
    958                          // odd
    959                          weights_[i] = 1.0;
    960                         //printf("bad pivot row %d (%d->%d) iPivot %d\n",i,firstRow,lastRow,iPivot);
    961                     }
    962                }
     896      double *array = savedWeights_->denseVector();
     897      int *which = savedWeights_->getIndices();
     898      for (int i = 0; i < numberRows; i++) {
     899        array[i] = weights_[i];
     900        which[i] = pivotVariable[i];
     901      }
     902      if (mode == 7) {
     903        savedWeights_->setNumElements(numberRows); //flag as special
     904        savedWeights_->setPackedMode(true);
     905      }
     906    } else if (mode != 6) {
     907      int *which = alternateWeights_->getIndices();
     908      CoinIndexedVector *rowArray3 = model_->rowArray(3);
     909      rowArray3->clear();
     910      int *back = rowArray3->getIndices();
     911      // In case something went wrong
     912      for (i = 0; i < numberRows + numberColumns; i++)
     913        back[i] = -1;
     914      if (mode != 4) {
     915        // save
     916        CoinMemcpyN(which, numberRows, savedWeights_->getIndices());
     917        CoinMemcpyN(weights_, numberRows, savedWeights_->denseVector());
     918      } else {
     919        // restore
     920        //memcpy(which,savedWeights_->getIndices(),
     921        //     numberRows*sizeof(int));
     922        //memcpy(weights_,savedWeights_->denseVector(),
     923        //     numberRows*sizeof(double));
     924        which = savedWeights_->getIndices();
     925      }
     926      // restore (a bit slow - but only every re-factorization)
     927      double *array = savedWeights_->denseVector();
     928      for (i = 0; i < numberRows; i++) {
     929        int iSeq = which[i];
     930        back[iSeq] = i;
     931      }
     932      int firstRow = 0;
     933      int lastRow = numberRows;
     934      if (mode == 7) {
     935        // use info passed in
     936        firstRow = model->spareIntArray_[0];
     937        lastRow = model->spareIntArray_[1];
     938      }
     939      for (i = firstRow; i < lastRow; i++) {
     940        int iPivot = pivotVariable[i];
     941        iPivot = back[iPivot];
     942        if (iPivot >= 0) {
     943          weights_[i] = array[iPivot];
     944          if (weights_[i] < DEVEX_TRY_NORM)
     945            weights_[i] = DEVEX_TRY_NORM; // may need to check more
     946        } else {
     947          // odd
     948          weights_[i] = 1.0;
     949          //printf("bad pivot row %d (%d->%d) iPivot %d\n",i,firstRow,lastRow,iPivot);
     950        }
     951      }
    963952#if 0
    964953               printf("mode %d mode_ %d state_ %d\n",mode,mode_,state_);
     
    969958               }
    970959#endif
    971           } else {
    972                // mode 6 - scale back weights as primal errors
    973                double primalError = model_->largestPrimalError();
    974                double allowed ;
    975                if (primalError > 1.0e3)
    976                     allowed = 10.0;
    977                else if (primalError > 1.0e2)
    978                     allowed = 50.0;
    979                else if (primalError > 1.0e1)
    980                     allowed = 100.0;
    981                else
    982                     allowed = 1000.0;
    983                double allowedInv = 1.0 / allowed;
    984                for (i = 0; i < numberRows; i++) {
    985                     double value = weights_[i];
    986                     if (value < allowedInv)
    987                          value = allowedInv;
    988                     else if (value > allowed)
    989                          value = allowed;
    990                     weights_[i] = allowed;
    991                }
    992           }
    993           state_ = 0;
    994           // set up infeasibilities
    995           if (!infeasible_) {
    996                infeasible_ = new CoinIndexedVector();
    997                infeasible_->reserve(numberRows);
    998           }
    999      }
    1000      if (mode >= 2) {
    1001           // Get dubious weights
    1002           //if (!dubiousWeights_)
    1003           //dubiousWeights_=new int[numberRows];
    1004           //model_->factorization()->getWeights(dubiousWeights_);
    1005           infeasible_->clear();
    1006           int iRow;
    1007           const int * pivotVariable = model_->pivotVariable();
    1008           double tolerance = model_->currentPrimalTolerance();
    1009           for (iRow = 0; iRow < numberRows; iRow++) {
    1010                int iPivot = pivotVariable[iRow];
    1011                double value = model_->solution(iPivot);
    1012                double lower = model_->lower(iPivot);
    1013                double upper = model_->upper(iPivot);
    1014                if (value < lower - tolerance) {
    1015                     value -= lower;
    1016                     value *= value;
     960    } else {
     961      // mode 6 - scale back weights as primal errors
     962      double primalError = model_->largestPrimalError();
     963      double allowed;
     964      if (primalError > 1.0e3)
     965        allowed = 10.0;
     966      else if (primalError > 1.0e2)
     967        allowed = 50.0;
     968      else if (primalError > 1.0e1)
     969        allowed = 100.0;
     970      else
     971        allowed = 1000.0;
     972      double allowedInv = 1.0 / allowed;
     973      for (i = 0; i < numberRows; i++) {
     974        double value = weights_[i];
     975        if (value < allowedInv)
     976          value = allowedInv;
     977        else if (value > allowed)
     978          value = allowed;
     979        weights_[i] = allowed;
     980      }
     981    }
     982    state_ = 0;
     983    // set up infeasibilities
     984    if (!infeasible_) {
     985      infeasible_ = new CoinIndexedVector();
     986      infeasible_->reserve(numberRows);
     987    }
     988  }
     989  if (mode >= 2) {
     990    // Get dubious weights
     991    //if (!dubiousWeights_)
     992    //dubiousWeights_=new int[numberRows];
     993    //model_->factorization()->getWeights(dubiousWeights_);
     994    infeasible_->clear();
     995    int iRow;
     996    const int *pivotVariable = model_->pivotVariable();
     997    double tolerance = model_->currentPrimalTolerance();
     998    for (iRow = 0; iRow < numberRows; iRow++) {
     999      int iPivot = pivotVariable[iRow];
     1000      double value = model_->solution(iPivot);
     1001      double lower = model_->lower(iPivot);
     1002      double upper = model_->upper(iPivot);
     1003      if (value < lower - tolerance) {
     1004        value -= lower;
     1005        value *= value;
    10171006#ifdef CLP_DUAL_COLUMN_MULTIPLIER
    1018                     if (iPivot < numberColumns)
    1019                          value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
     1007        if (iPivot < numberColumns)
     1008          value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
    10201009#endif
    10211010#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
    1022                     if (lower == upper)
    1023                          value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
    1024 #endif
    1025                     // store square in list
    1026                     infeasible_->quickAdd(iRow, value);
    1027                } else if (value > upper + tolerance) {
    1028                     value -= upper;
    1029                     value *= value;
     1011        if (lower == upper)
     1012          value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
     1013#endif
     1014        // store square in list
     1015        infeasible_->quickAdd(iRow, value);
     1016      } else if (value > upper + tolerance) {
     1017        value -= upper;
     1018        value *= value;
    10301019#ifdef CLP_DUAL_COLUMN_MULTIPLIER
    1031                     if (iPivot < numberColumns)
    1032                          value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
     1020        if (iPivot < numberColumns)
     1021          value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
    10331022#endif
    10341023#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
    1035                     if (lower == upper)
    1036                          value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
    1037 #endif
    1038                     // store square in list
    1039                     infeasible_->quickAdd(iRow, value);
    1040                }
    1041           }
    1042      }
    1043      // see where coming from
    1044      if (mode==2&&!model->numberIterations()) {
    1045        int options=model->specialOptions();
    1046        if ((options&16384)!=0&&true) {
    1047         // fast of some sort - could be clever???
    1048         // for now initialize
    1049          if ((options&524288)!=0&&false) {
    1050            // fathom
    1051            for (int i = 0; i < numberRows; i++)
    1052              weights_[i] = 1.0;
    1053         } else if (true) {
    1054            // strong
    1055            for (int i = 0; i < numberRows; i++)
    1056              weights_[i] = 1.0;
    1057         }
    1058        }
    1059      }
     1024        if (lower == upper)
     1025          value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
     1026#endif
     1027        // store square in list
     1028        infeasible_->quickAdd(iRow, value);
     1029      }
     1030    }
     1031  }
     1032  // see where coming from
     1033  if (mode == 2 && !model->numberIterations()) {
     1034    int options = model->specialOptions();
     1035    if ((options & 16384) != 0 && true) {
     1036      // fast of some sort - could be clever???
     1037      // for now initialize
     1038      if ((options & 524288) != 0 && false) {
     1039        // fathom
     1040        for (int i = 0; i < numberRows; i++)
     1041          weights_[i] = 1.0;
     1042      } else if (true) {
     1043        // strong
     1044        for (int i = 0; i < numberRows; i++)
     1045          weights_[i] = 1.0;
     1046      }
     1047    }
     1048  }
    10601049}
    10611050// Pass in saved weights
    1062 void
    1063 ClpDualRowSteepest::passInSavedWeights(const CoinIndexedVector * saved)
     1051void ClpDualRowSteepest::passInSavedWeights(const CoinIndexedVector *saved)
    10641052{
    10651053  delete savedWeights_;
    1066   savedWeights_=new CoinIndexedVector(*saved);
     1054  savedWeights_ = new CoinIndexedVector(*saved);
    10671055}
    10681056// Gets rid of last update
    1069 void
    1070 ClpDualRowSteepest::unrollWeights()
    1071 {
    1072      double * saved = alternateWeights_->denseVector();
    1073      int number = alternateWeights_->getNumElements();
    1074      int * which = alternateWeights_->getIndices();
    1075      int i;
    1076      if (alternateWeights_->packedMode()) {
    1077           for (i = 0; i < number; i++) {
    1078                int iRow = which[i];
    1079                weights_[iRow] = saved[i];
    1080                saved[i] = 0.0;
    1081           }
    1082      } else {
    1083           for (i = 0; i < number; i++) {
    1084                int iRow = which[i];
    1085                weights_[iRow] = saved[iRow];
    1086                saved[iRow] = 0.0;
    1087           }
    1088      }
    1089      alternateWeights_->setNumElements(0);
     1057void ClpDualRowSteepest::unrollWeights()
     1058{
     1059  double *saved = alternateWeights_->denseVector();
     1060  int number = alternateWeights_->getNumElements();
     1061  int *which = alternateWeights_->getIndices();
     1062  int i;
     1063  if (alternateWeights_->packedMode()) {
     1064    for (i = 0; i < number; i++) {
     1065      int iRow = which[i];
     1066      weights_[iRow] = saved[i];
     1067      saved[i] = 0.0;
     1068    }
     1069  } else {
     1070    for (i = 0; i < number; i++) {
     1071      int iRow = which[i];
     1072      weights_[iRow] = saved[iRow];
     1073      saved[iRow] = 0.0;
     1074    }
     1075  }
     1076  alternateWeights_->setNumElements(0);
    10901077}
    10911078//-------------------------------------------------------------------
    10921079// Clone
    10931080//-------------------------------------------------------------------
    1094 ClpDualRowPivot * ClpDualRowSteepest::clone(bool CopyData) const
    1095 {
    1096      if (CopyData) {
    1097           return new ClpDualRowSteepest(*this);
    1098      } else {
    1099           return new ClpDualRowSteepest();
    1100      }
     1081ClpDualRowPivot *ClpDualRowSteepest::clone(bool CopyData) const
     1082{
     1083  if (CopyData) {
     1084    return new ClpDualRowSteepest(*this);
     1085  } else {
     1086    return new ClpDualRowSteepest();
     1087  }
    11011088}
    11021089// Gets rid of all arrays
    1103 void
    1104 ClpDualRowSteepest::clearArrays()
    1105 {
    1106      if (persistence_ == normal) {
    1107           delete [] weights_;
    1108           weights_ = NULL;
    1109           delete [] dubiousWeights_;
    1110           dubiousWeights_ = NULL;
    1111           delete infeasible_;
    1112           infeasible_ = NULL;
    1113           delete alternateWeights_;
    1114           alternateWeights_ = NULL;
    1115           delete savedWeights_;
    1116           savedWeights_ = NULL;
    1117      }
    1118      state_ = -1;
     1090void ClpDualRowSteepest::clearArrays()
     1091{
     1092  if (persistence_ == normal) {
     1093    delete[] weights_;
     1094    weights_ = NULL;
     1095    delete[] dubiousWeights_;
     1096    dubiousWeights_ = NULL;
     1097    delete infeasible_;
     1098    infeasible_ = NULL;
     1099    delete alternateWeights_;
     1100    alternateWeights_ = NULL;
     1101    delete savedWeights_;
     1102    savedWeights_ = NULL;
     1103  }
     1104  state_ = -1;
    11191105}
    11201106// Returns true if would not find any row
    1121 bool
    1122 ClpDualRowSteepest::looksOptimal() const
    1123 {
    1124      int iRow;
    1125      const int * pivotVariable = model_->pivotVariable();
    1126      double tolerance = model_->currentPrimalTolerance();
    1127      // we can't really trust infeasibilities if there is primal error
    1128      // this coding has to mimic coding in checkPrimalSolution
    1129      double error = CoinMin(1.0e-2, model_->largestPrimalError());
    1130      // allow tolerance at least slightly bigger than standard
    1131      tolerance = tolerance +  error;
    1132      // But cap
    1133      tolerance = CoinMin(1000.0, tolerance);
    1134      int numberRows = model_->numberRows();
    1135      int numberInfeasible = 0;
    1136      for (iRow = 0; iRow < numberRows; iRow++) {
    1137           int iPivot = pivotVariable[iRow];
    1138           double value = model_->solution(iPivot);
    1139           double lower = model_->lower(iPivot);
    1140           double upper = model_->upper(iPivot);
    1141           if (value < lower - tolerance) {
    1142                numberInfeasible++;
    1143           } else if (value > upper + tolerance) {
    1144                numberInfeasible++;
    1145           }
    1146      }
    1147      return (numberInfeasible == 0);
     1107bool ClpDualRowSteepest::looksOptimal() const
     1108{
     1109  int iRow;
     1110  const int *pivotVariable = model_->pivotVariable();
     1111  double tolerance = model_->currentPrimalTolerance();
     1112  // we can't really trust infeasibilities if there is primal error
     1113  // this coding has to mimic coding in checkPrimalSolution
     1114  double error = CoinMin(1.0e-2, model_->largestPrimalError());
     1115  // allow tolerance at least slightly bigger than standard
     1116  tolerance = tolerance + error;
     1117  // But cap
     1118  tolerance = CoinMin(1000.0, tolerance);
     1119  int numberRows = model_->numberRows();
     1120  int numberInfeasible = 0;
     1121  for (iRow = 0; iRow < numberRows; iRow++) {
     1122    int iPivot = pivotVariable[iRow];
     1123    double value = model_->solution(iPivot);
     1124    double lower = model_->lower(iPivot);
     1125    double upper = model_->upper(iPivot);
     1126    if (value < lower - tolerance) {
     1127      numberInfeasible++;
     1128    } else if (value > upper + tolerance) {
     1129      numberInfeasible++;
     1130    }
     1131  }
     1132  return (numberInfeasible == 0);
    11481133}
    11491134// Called when maximum pivots changes
    1150 void
    1151 ClpDualRowSteepest::maximumPivotsChanged()
    1152 {
    1153      if (alternateWeights_ &&
    1154                alternateWeights_->capacity() != model_->numberRows() +
    1155                model_->factorization()->maximumPivots()) {
    1156           delete alternateWeights_;
    1157           alternateWeights_ = new CoinIndexedVector();
    1158           // enough space so can use it for factorization
    1159           alternateWeights_->reserve(model_->numberRows() +
    1160                                      model_->factorization()->maximumPivots());
    1161      }
     1135void ClpDualRowSteepest::maximumPivotsChanged()
     1136{
     1137  if (alternateWeights_ && alternateWeights_->capacity() != model_->numberRows() + model_->factorization()->maximumPivots()) {
     1138    delete alternateWeights_;
     1139    alternateWeights_ = new CoinIndexedVector();
     1140    // enough space so can use it for factorization
     1141    alternateWeights_->reserve(model_->numberRows() + model_->factorization()->maximumPivots());
     1142  }
    11621143}
    11631144
     1145/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
     1146*/
Note: See TracChangeset for help on using the changeset viewer.