/* $Id: ClpDualRowSteepest.cpp 1732 2011-05-31 08:09:41Z forrest $ */
// Copyright (C) 2002, International Business Machines
// Corporation and others.  All Rights Reserved.
// This code is licensed under the terms of the Eclipse Public License (EPL).

#include "CoinPragma.hpp"
#include "ClpSimplex.hpp"
#include "ClpDualRowSteepest.hpp"
#include "CoinIndexedVector.hpp"
#include "ClpFactorization.hpp"
#include "CoinHelperFunctions.hpp"
#include <cstdio>
//#############################################################################
// Constructors / Destructor / Assignment
//#############################################################################
//#define CLP_DEBUG 4
//-------------------------------------------------------------------
// Default Constructor
//-------------------------------------------------------------------
ClpDualRowSteepest::ClpDualRowSteepest (int mode)
     : ClpDualRowPivot(),
       state_(-1),
       mode_(mode),
       persistence_(normal),
       weights_(NULL),
       infeasible_(NULL),
       alternateWeights_(NULL),
       savedWeights_(NULL),
       dubiousWeights_(NULL)
{
     type_ = 2 + 64 * mode;
}

//-------------------------------------------------------------------
// Copy constructor
//-------------------------------------------------------------------
ClpDualRowSteepest::ClpDualRowSteepest (const ClpDualRowSteepest & rhs)
     : ClpDualRowPivot(rhs)
{
     state_ = rhs.state_;
     mode_ = rhs.mode_;
     persistence_ = rhs.persistence_;
     model_ = rhs.model_;
     if ((model_ && model_->whatsChanged() & 1) != 0) {
          int number = model_->numberRows();
          if (rhs.savedWeights_)
               number = CoinMin(number, rhs.savedWeights_->capacity());
          if (rhs.infeasible_) {
               infeasible_ = new CoinIndexedVector(rhs.infeasible_);
          } else {
               infeasible_ = NULL;
          }
          if (rhs.weights_) {
               weights_ = new double[number];
               ClpDisjointCopyN(rhs.weights_, number, weights_);
          } else {
               weights_ = NULL;
          }
          if (rhs.alternateWeights_) {
               alternateWeights_ = new CoinIndexedVector(rhs.alternateWeights_);
          } else {
               alternateWeights_ = NULL;
          }
          if (rhs.savedWeights_) {
               savedWeights_ = new CoinIndexedVector(rhs.savedWeights_);
          } else {
               savedWeights_ = NULL;
          }
          if (rhs.dubiousWeights_) {
               assert(model_);
               int number = model_->numberRows();
               dubiousWeights_ = new int[number];
               ClpDisjointCopyN(rhs.dubiousWeights_, number, dubiousWeights_);
          } else {
               dubiousWeights_ = NULL;
          }
     } else {
          infeasible_ = NULL;
          weights_ = NULL;
          alternateWeights_ = NULL;
          savedWeights_ = NULL;
          dubiousWeights_ = NULL;
     }
}

//-------------------------------------------------------------------
// Destructor
//-------------------------------------------------------------------
ClpDualRowSteepest::~ClpDualRowSteepest ()
{
     delete [] weights_;
     delete [] dubiousWeights_;
     delete infeasible_;
     delete alternateWeights_;
     delete savedWeights_;
}

//----------------------------------------------------------------
// Assignment operator
//-------------------------------------------------------------------
ClpDualRowSteepest &
ClpDualRowSteepest::operator=(const ClpDualRowSteepest& rhs)
{
     if (this != &rhs) {
          ClpDualRowPivot::operator=(rhs);
          state_ = rhs.state_;
          mode_ = rhs.mode_;
          persistence_ = rhs.persistence_;
          model_ = rhs.model_;
          delete [] weights_;
          delete [] dubiousWeights_;
          delete infeasible_;
          delete alternateWeights_;
          delete savedWeights_;
          assert(model_);
          int number = model_->numberRows();
          if (rhs.savedWeights_)
               number = CoinMin(number, rhs.savedWeights_->capacity());
          if (rhs.infeasible_ != NULL) {
               infeasible_ = new CoinIndexedVector(rhs.infeasible_);
          } else {
               infeasible_ = NULL;
          }
          if (rhs.weights_ != NULL) {
               weights_ = new double[number];
               ClpDisjointCopyN(rhs.weights_, number, weights_);
          } else {
               weights_ = NULL;
          }
          if (rhs.alternateWeights_ != NULL) {
               alternateWeights_ = new CoinIndexedVector(rhs.alternateWeights_);
          } else {
               alternateWeights_ = NULL;
          }
          if (rhs.savedWeights_ != NULL) {
               savedWeights_ = new CoinIndexedVector(rhs.savedWeights_);
          } else {
               savedWeights_ = NULL;
          }
          if (rhs.dubiousWeights_) {
               assert(model_);
               int number = model_->numberRows();
               dubiousWeights_ = new int[number];
               ClpDisjointCopyN(rhs.dubiousWeights_, number, dubiousWeights_);
          } else {
               dubiousWeights_ = NULL;
          }
     }
     return *this;
}
// Fill most values
void
ClpDualRowSteepest::fill(const ClpDualRowSteepest& rhs)
{
     state_ = rhs.state_;
     mode_ = rhs.mode_;
     persistence_ = rhs.persistence_;
     assert (model_->numberRows() == rhs.model_->numberRows());
     model_ = rhs.model_;
     assert(model_);
     int number = model_->numberRows();
     if (rhs.savedWeights_)
          number = CoinMin(number, rhs.savedWeights_->capacity());
     if (rhs.infeasible_ != NULL) {
          if (!infeasible_)
               infeasible_ = new CoinIndexedVector(rhs.infeasible_);
          else
               *infeasible_ = *rhs.infeasible_;
     } else {
          delete infeasible_;
          infeasible_ = NULL;
     }
     if (rhs.weights_ != NULL) {
          if (!weights_)
               weights_ = new double[number];
          ClpDisjointCopyN(rhs.weights_, number, weights_);
     } else {
          delete [] weights_;
          weights_ = NULL;
     }
     if (rhs.alternateWeights_ != NULL) {
          if (!alternateWeights_)
               alternateWeights_ = new CoinIndexedVector(rhs.alternateWeights_);
          else
               *alternateWeights_ = *rhs.alternateWeights_;
     } else {
          delete alternateWeights_;
          alternateWeights_ = NULL;
     }
     if (rhs.savedWeights_ != NULL) {
          if (!savedWeights_)
               savedWeights_ = new CoinIndexedVector(rhs.savedWeights_);
          else
               *savedWeights_ = *rhs.savedWeights_;
     } else {
          delete savedWeights_;
          savedWeights_ = NULL;
     }
     if (rhs.dubiousWeights_) {
          assert(model_);
          int number = model_->numberRows();
          if (!dubiousWeights_)
               dubiousWeights_ = new int[number];
          ClpDisjointCopyN(rhs.dubiousWeights_, number, dubiousWeights_);
     } else {
          delete [] dubiousWeights_;
          dubiousWeights_ = NULL;
     }
}
// Returns pivot row, -1 if none
int
ClpDualRowSteepest::pivotRow()
{
     assert(model_);
     int i, iRow;
     double * infeas = infeasible_->denseVector();
     double largest = 0.0;
     int * index = infeasible_->getIndices();
     int number = infeasible_->getNumElements();
     const int * pivotVariable = model_->pivotVariable();
     int chosenRow = -1;
     int lastPivotRow = model_->pivotRow();
     assert (lastPivotRow < model_->numberRows());
     double tolerance = model_->currentPrimalTolerance();
     // we can't really trust infeasibilities if there is primal error
     // this coding has to mimic coding in checkPrimalSolution
     double error = CoinMin(1.0e-2, model_->largestPrimalError());
     // allow tolerance at least slightly bigger than standard
     tolerance = tolerance +  error;
     // But cap
     tolerance = CoinMin(1000.0, tolerance);
     tolerance *= tolerance; // as we are using squares
     bool toleranceChanged = false;
     double * solution = model_->solutionRegion();
     double * lower = model_->lowerRegion();
     double * upper = model_->upperRegion();
     // do last pivot row one here
     //#define CLP_DUAL_FIXED_COLUMN_MULTIPLIER 10.0
     if (lastPivotRow >= 0 && lastPivotRow < model_->numberRows()) {
#ifdef CLP_DUAL_COLUMN_MULTIPLIER
          int numberColumns = model_->numberColumns();
#endif
          int iPivot = pivotVariable[lastPivotRow];
          double value = solution[iPivot];
          double lower = model_->lower(iPivot);
          double upper = model_->upper(iPivot);
          if (value > upper + tolerance) {
               value -= upper;
               value *= value;
#ifdef CLP_DUAL_COLUMN_MULTIPLIER
               if (iPivot < numberColumns)
                    value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
#endif
               // store square in list
               if (infeas[lastPivotRow])
                    infeas[lastPivotRow] = value; // already there
               else
                    infeasible_->quickAdd(lastPivotRow, value);
          } else if (value < lower - tolerance) {
               value -= lower;
               value *= value;
#ifdef CLP_DUAL_COLUMN_MULTIPLIER
               if (iPivot < numberColumns)
                    value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
#endif
               // store square in list
               if (infeas[lastPivotRow])
                    infeas[lastPivotRow] = value; // already there
               else
                    infeasible_->add(lastPivotRow, value);
          } else {
               // feasible - was it infeasible - if so set tiny
               if (infeas[lastPivotRow])
                    infeas[lastPivotRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
          }
          number = infeasible_->getNumElements();
     }
     if(model_->numberIterations() < model_->lastBadIteration() + 200) {
          // we can't really trust infeasibilities if there is dual error
          if (model_->largestDualError() > model_->largestPrimalError()) {
               tolerance *= CoinMin(model_->largestDualError() / model_->largestPrimalError(), 1000.0);
               toleranceChanged = true;
          }
     }
     int numberWanted;
     if (mode_ < 2 ) {
          numberWanted = number + 1;
     } else if (mode_ == 2) {
          numberWanted = CoinMax(2000, number / 8);
     } else {
          int numberElements = model_->factorization()->numberElements();
          double ratio = static_cast<double> (numberElements) /
                         static_cast<double> (model_->numberRows());
          numberWanted = CoinMax(2000, number / 8);
          if (ratio < 1.0) {
               numberWanted = CoinMax(2000, number / 20);
          } else if (ratio > 10.0) {
               ratio = number * (ratio / 80.0);
               if (ratio > number)
                    numberWanted = number + 1;
               else
                    numberWanted = CoinMax(2000, static_cast<int> (ratio));
          }
     }
     if (model_->largestPrimalError() > 1.0e-3)
          numberWanted = number + 1; // be safe
     int iPass;
     // Setup two passes
     int start[4];
     start[1] = number;
     start[2] = 0;
     double dstart = static_cast<double> (number) * model_->randomNumberGenerator()->randomDouble();
     start[0] = static_cast<int> (dstart);
     start[3] = start[0];
     //double largestWeight=0.0;
     //double smallestWeight=1.0e100;
     for (iPass = 0; iPass < 2; iPass++) {
          int end = start[2*iPass+1];
          for (i = start[2*iPass]; i < end; i++) {
               iRow = index[i];
               double value = infeas[iRow];
               if (value > tolerance) {
                    //#define OUT_EQ
#ifdef OUT_EQ
                    {
                         int iSequence = pivotVariable[iRow];
                         if (upper[iSequence] == lower[iSequence])
                              value *= 2.0;
                    }
#endif
                    double weight = CoinMin(weights_[iRow], 1.0e50);
                    //largestWeight = CoinMax(largestWeight,weight);
                    //smallestWeight = CoinMin(smallestWeight,weight);
                    //double dubious = dubiousWeights_[iRow];
                    //weight *= dubious;
                    //if (value>2.0*largest*weight||(value>0.5*largest*weight&&value*largestWeight>dubious*largest*weight)) {
                    if (value > largest * weight) {
                         // make last pivot row last resort choice
                         if (iRow == lastPivotRow) {
                              if (value * 1.0e-10 < largest * weight)
                                   continue;
                              else
                                   value *= 1.0e-10;
                         }
                         int iSequence = pivotVariable[iRow];
                         if (!model_->flagged(iSequence)) {
                              //#define CLP_DEBUG 3
#ifdef CLP_DEBUG
                              double value2 = 0.0;
                              if (solution[iSequence] > upper[iSequence] + tolerance)
                                   value2 = solution[iSequence] - upper[iSequence];
                              else if (solution[iSequence] < lower[iSequence] - tolerance)
                                   value2 = solution[iSequence] - lower[iSequence];
                              assert(fabs(value2 * value2 - infeas[iRow]) < 1.0e-8 * CoinMin(value2 * value2, infeas[iRow]));
#endif
                              if (solution[iSequence] > upper[iSequence] + tolerance ||
                                        solution[iSequence] < lower[iSequence] - tolerance) {
                                   chosenRow = iRow;
                                   largest = value / weight;
                                   //largestWeight = dubious;
                              }
                         } else {
                              // just to make sure we don't exit before got something
                              numberWanted++;
                         }
                    }
                    numberWanted--;
                    if (!numberWanted)
                         break;
               }
          }
          if (!numberWanted)
               break;
     }
     //printf("smallest %g largest %g\n",smallestWeight,largestWeight);
     if (chosenRow < 0 && toleranceChanged) {
          // won't line up with checkPrimalSolution - do again
          double saveError = model_->largestDualError();
          model_->setLargestDualError(0.0);
          // can't loop
          chosenRow = pivotRow();
          model_->setLargestDualError(saveError);
     }
     if (chosenRow < 0 && lastPivotRow < 0) {
          int nLeft = 0;
          for (int i = 0; i < number; i++) {
               int iRow = index[i];
               if (fabs(infeas[iRow]) > 1.0e-50) {
                    index[nLeft++] = iRow;
               } else {
                    infeas[iRow] = 0.0;
               }
          }
          infeasible_->setNumElements(nLeft);
          model_->setNumberPrimalInfeasibilities(nLeft);
     }
     return chosenRow;
}
#if 0
static double ft_count = 0.0;
static double up_count = 0.0;
static double ft_count_in = 0.0;
static double up_count_in = 0.0;
static int xx_count = 0;
#endif
/* Updates weights and returns pivot alpha.
   Also does FT update */
double
ClpDualRowSteepest::updateWeights(CoinIndexedVector * input,
                                  CoinIndexedVector * spare,
                                  CoinIndexedVector * spare2,
                                  CoinIndexedVector * updatedColumn)
{
     //#define CLP_DEBUG 3
#if CLP_DEBUG>2
     // Very expensive debug
     {
          int numberRows = model_->numberRows();
          CoinIndexedVector * temp = new CoinIndexedVector();
          temp->reserve(numberRows +
                        model_->factorization()->maximumPivots());
          double * array = alternateWeights_->denseVector();
          int * which = alternateWeights_->getIndices();
          alternateWeights_->clear();
          int i;
          for (i = 0; i < numberRows; i++) {
               double value = 0.0;
               array[i] = 1.0;
               which[0] = i;
               alternateWeights_->setNumElements(1);
               model_->factorization()->updateColumnTranspose(temp,
                         alternateWeights_);
               int number = alternateWeights_->getNumElements();
               int j;
               for (j = 0; j < number; j++) {
                    int iRow = which[j];
                    value += array[iRow] * array[iRow];
                    array[iRow] = 0.0;
               }
               alternateWeights_->setNumElements(0);
               double w = CoinMax(weights_[i], value) * .1;
               if (fabs(weights_[i] - value) > w) {
                    printf("%d old %g, true %g\n", i, weights_[i], value);
                    weights_[i] = value; // to reduce printout
               }
               //else
               //printf("%d matches %g\n",i,value);
          }
          delete temp;
     }
#endif
     assert (input->packedMode());
     if (!updatedColumn->packedMode()) {
          // I think this means empty
#ifdef COIN_DEVELOP
          printf("updatedColumn not packed mode ClpDualRowSteepest::updateWeights\n");
#endif
          return 0.0;
     }
     double alpha = 0.0;
     if (!model_->factorization()->networkBasis()) {
          // clear other region
          alternateWeights_->clear();
          double norm = 0.0;
          int i;
          double * work = input->denseVector();
          int numberNonZero = input->getNumElements();
          int * which = input->getIndices();
          double * work2 = spare->denseVector();
          int * which2 = spare->getIndices();
          // ftran
          //permute and move indices into index array
          //also compute norm
          //int *regionIndex = alternateWeights_->getIndices (  );
          const int *permute = model_->factorization()->permute();
          //double * region = alternateWeights_->denseVector();
          if (permute) {
               for ( i = 0; i < numberNonZero; i ++ ) {
                    int iRow = which[i];
                    double value = work[i];
                    norm += value * value;
                    iRow = permute[iRow];
                    work2[iRow] = value;
                    which2[i] = iRow;
               }
          } else {
               for ( i = 0; i < numberNonZero; i ++ ) {
                    int iRow = which[i];
                    double value = work[i];
                    norm += value * value;
                    //iRow = permute[iRow];
                    work2[iRow] = value;
                    which2[i] = iRow;
               }
          }
          spare->setNumElements ( numberNonZero );
          // Do FT update
#if 0
          ft_count_in += updatedColumn->getNumElements();
          up_count_in += spare->getNumElements();
#endif
          if (permute || true) {
#if CLP_DEBUG>2
               printf("REGION before %d els\n", spare->getNumElements());
               spare->print();
#endif
               model_->factorization()->updateTwoColumnsFT(spare2, updatedColumn,
                         spare, permute != NULL);
#if CLP_DEBUG>2
               printf("REGION after %d els\n", spare->getNumElements());
               spare->print();
#endif
          } else {
               // Leave as old way
               model_->factorization()->updateColumnFT(spare2, updatedColumn);
               model_->factorization()->updateColumn(spare2, spare, false);
          }
#undef CLP_DEBUG
#if 0
          ft_count += updatedColumn->getNumElements();
          up_count += spare->getNumElements();
          xx_count++;
          if ((xx_count % 1000) == 0)
               printf("zz %d ft %g up %g (in %g %g)\n", xx_count, ft_count, up_count,
                      ft_count_in, up_count_in);
#endif
          numberNonZero = spare->getNumElements();
          // alternateWeights_ should still be empty
          int pivotRow = model_->pivotRow();
#ifdef CLP_DEBUG
          if ( model_->logLevel (  ) > 4  &&
                    fabs(norm - weights_[pivotRow]) > 1.0e-3 * (1.0 + norm))
               printf("on row %d, true weight %g, old %g\n",
                      pivotRow, sqrt(norm), sqrt(weights_[pivotRow]));
#endif
          // could re-initialize here (could be expensive)
          norm /= model_->alpha() * model_->alpha();
          assert(model_->alpha());
          assert(norm);
          // pivot element
          alpha = 0.0;
          double multiplier = 2.0 / model_->alpha();
          // look at updated column
          work = updatedColumn->denseVector();
          numberNonZero = updatedColumn->getNumElements();
          which = updatedColumn->getIndices();

          int nSave = 0;
          double * work3 = alternateWeights_->denseVector();
          int * which3 = alternateWeights_->getIndices();
          const int * pivotColumn = model_->factorization()->pivotColumn();
          for (i = 0; i < numberNonZero; i++) {
               int iRow = which[i];
               double theta = work[i];
               if (iRow == pivotRow)
                    alpha = theta;
               double devex = weights_[iRow];
               work3[nSave] = devex; // save old
               which3[nSave++] = iRow;
               // transform to match spare
               int jRow = permute ? pivotColumn[iRow] : iRow;
               double value = work2[jRow];
               devex +=  theta * (theta * norm + value * multiplier);
               if (devex < DEVEX_TRY_NORM)
                    devex = DEVEX_TRY_NORM;
               weights_[iRow] = devex;
          }
          alternateWeights_->setPackedMode(true);
          alternateWeights_->setNumElements(nSave);
          if (norm < DEVEX_TRY_NORM)
               norm = DEVEX_TRY_NORM;
          // Try this to make less likely will happen again and stop cycling
          //norm *= 1.02;
          weights_[pivotRow] = norm;
          spare->clear();
#ifdef CLP_DEBUG
          spare->checkClear();
#endif
     } else {
          // Do FT update
          model_->factorization()->updateColumnFT(spare, updatedColumn);
          // clear other region
          alternateWeights_->clear();
          double norm = 0.0;
          int i;
          double * work = input->denseVector();
          int number = input->getNumElements();
          int * which = input->getIndices();
          double * work2 = spare->denseVector();
          int * which2 = spare->getIndices();
          for (i = 0; i < number; i++) {
               int iRow = which[i];
               double value = work[i];
               norm += value * value;
               work2[iRow] = value;
               which2[i] = iRow;
          }
          spare->setNumElements(number);
          // ftran
#ifndef NDEBUG
          alternateWeights_->checkClear();
#endif
          model_->factorization()->updateColumn(alternateWeights_, spare);
          // alternateWeights_ should still be empty
#ifndef NDEBUG
          alternateWeights_->checkClear();
#endif
          int pivotRow = model_->pivotRow();
#ifdef CLP_DEBUG
          if ( model_->logLevel (  ) > 4  &&
                    fabs(norm - weights_[pivotRow]) > 1.0e-3 * (1.0 + norm))
               printf("on row %d, true weight %g, old %g\n",
                      pivotRow, sqrt(norm), sqrt(weights_[pivotRow]));
#endif
          // could re-initialize here (could be expensive)
          norm /= model_->alpha() * model_->alpha();

          assert(norm);
          //if (norm < DEVEX_TRY_NORM)
          //norm = DEVEX_TRY_NORM;
          // pivot element
          alpha = 0.0;
          double multiplier = 2.0 / model_->alpha();
          // look at updated column
          work = updatedColumn->denseVector();
          number = updatedColumn->getNumElements();
          which = updatedColumn->getIndices();

          int nSave = 0;
          double * work3 = alternateWeights_->denseVector();
          int * which3 = alternateWeights_->getIndices();
          for (i = 0; i < number; i++) {
               int iRow = which[i];
               double theta = work[i];
               if (iRow == pivotRow)
                    alpha = theta;
               double devex = weights_[iRow];
               work3[nSave] = devex; // save old
               which3[nSave++] = iRow;
               double value = work2[iRow];
               devex +=  theta * (theta * norm + value * multiplier);
               if (devex < DEVEX_TRY_NORM)
                    devex = DEVEX_TRY_NORM;
               weights_[iRow] = devex;
          }
          if (!alpha) {
               // error - but carry on
               alpha = 1.0e-50;
          }
          alternateWeights_->setPackedMode(true);
          alternateWeights_->setNumElements(nSave);
          if (norm < DEVEX_TRY_NORM)
               norm = DEVEX_TRY_NORM;
          weights_[pivotRow] = norm;
          spare->clear();
     }
#ifdef CLP_DEBUG
     spare->checkClear();
#endif
     return alpha;
}

/* Updates primal solution (and maybe list of candidates)
   Uses input vector which it deletes
   Computes change in objective function
*/
void
ClpDualRowSteepest::updatePrimalSolution(
     CoinIndexedVector * primalUpdate,
     double primalRatio,
     double & objectiveChange)
{
     double * COIN_RESTRICT work = primalUpdate->denseVector();
     int number = primalUpdate->getNumElements();
     int * COIN_RESTRICT which = primalUpdate->getIndices();
     int i;
     double changeObj = 0.0;
     double tolerance = model_->currentPrimalTolerance();
     const int * COIN_RESTRICT pivotVariable = model_->pivotVariable();
     double * COIN_RESTRICT infeas = infeasible_->denseVector();
     double * COIN_RESTRICT solution = model_->solutionRegion();
     const double * COIN_RESTRICT costModel = model_->costRegion();
     const double * COIN_RESTRICT lowerModel = model_->lowerRegion();
     const double * COIN_RESTRICT upperModel = model_->upperRegion();
#ifdef CLP_DUAL_COLUMN_MULTIPLIER
     int numberColumns = model_->numberColumns();
#endif
     if (primalUpdate->packedMode()) {
          for (i = 0; i < number; i++) {
               int iRow = which[i];
               int iPivot = pivotVariable[iRow];
               double value = solution[iPivot];
               double cost = costModel[iPivot];
               double change = primalRatio * work[i];
               work[i] = 0.0;
               value -= change;
               changeObj -= change * cost;
               double lower = lowerModel[iPivot];
               double upper = upperModel[iPivot];
               solution[iPivot] = value;
               if (value < lower - tolerance) {
                    value -= lower;
                    value *= value;
#ifdef CLP_DUAL_COLUMN_MULTIPLIER
                    if (iPivot < numberColumns)
                         value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
#endif
#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
                    if (lower == upper)
                         value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
#endif
                    // store square in list
                    if (infeas[iRow])
                         infeas[iRow] = value; // already there
                    else
                         infeasible_->quickAdd(iRow, value);
               } else if (value > upper + tolerance) {
                    value -= upper;
                    value *= value;
#ifdef CLP_DUAL_COLUMN_MULTIPLIER
                    if (iPivot < numberColumns)
                         value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
#endif
#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
                    if (lower == upper)
                         value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
#endif
                    // store square in list
                    if (infeas[iRow])
                         infeas[iRow] = value; // already there
                    else
                         infeasible_->quickAdd(iRow, value);
               } else {
                    // feasible - was it infeasible - if so set tiny
                    if (infeas[iRow])
                         infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
               }
          }
     } else {
          for (i = 0; i < number; i++) {
               int iRow = which[i];
               int iPivot = pivotVariable[iRow];
               double value = solution[iPivot];
               double cost = costModel[iPivot];
               double change = primalRatio * work[iRow];
               value -= change;
               changeObj -= change * cost;
               double lower = lowerModel[iPivot];
               double upper = upperModel[iPivot];
               solution[iPivot] = value;
               if (value < lower - tolerance) {
                    value -= lower;
                    value *= value;
#ifdef CLP_DUAL_COLUMN_MULTIPLIER
                    if (iPivot < numberColumns)
                         value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
#endif
#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
                    if (lower == upper)
                         value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
#endif
                    // store square in list
                    if (infeas[iRow])
                         infeas[iRow] = value; // already there
                    else
                         infeasible_->quickAdd(iRow, value);
               } else if (value > upper + tolerance) {
                    value -= upper;
                    value *= value;
#ifdef CLP_DUAL_COLUMN_MULTIPLIER
                    if (iPivot < numberColumns)
                         value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
#endif
#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
                    if (lower == upper)
                         value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
#endif
                    // store square in list
                    if (infeas[iRow])
                         infeas[iRow] = value; // already there
                    else
                         infeasible_->quickAdd(iRow, value);
               } else {
                    // feasible - was it infeasible - if so set tiny
                    if (infeas[iRow])
                         infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
               }
               work[iRow] = 0.0;
          }
     }
     // Do pivot row
     {
          int iRow = model_->pivotRow();
          // feasible - was it infeasible - if so set tiny
          //assert (infeas[iRow]);
          if (infeas[iRow])
               infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
     }
     primalUpdate->setNumElements(0);
     objectiveChange += changeObj;
}
/* Saves any weights round factorization as pivot rows may change
   1) before factorization
   2) after factorization
   3) just redo infeasibilities
   4) restore weights
*/
void
ClpDualRowSteepest::saveWeights(ClpSimplex * model, int mode)
{
     // alternateWeights_ is defined as indexed but is treated oddly
     model_ = model;
     int numberRows = model_->numberRows();
     int numberColumns = model_->numberColumns();
     const int * pivotVariable = model_->pivotVariable();
     int i;
     if (mode == 1) {
          if(weights_) {
               // Check if size has changed
               if (infeasible_->capacity() == numberRows) {
                    alternateWeights_->clear();
                    // change from row numbers to sequence numbers
                    int * which = alternateWeights_->getIndices();
                    for (i = 0; i < numberRows; i++) {
                         int iPivot = pivotVariable[i];
                         which[i] = iPivot;
                    }
                    state_ = 1;
               } else {
                    // size has changed - clear everything
                    delete [] weights_;
                    weights_ = NULL;
                    delete [] dubiousWeights_;
                    dubiousWeights_ = NULL;
                    delete infeasible_;
                    infeasible_ = NULL;
                    delete alternateWeights_;
                    alternateWeights_ = NULL;
                    delete savedWeights_;
                    savedWeights_ = NULL;
                    state_ = -1;
               }
          }
     } else if (mode == 2 || mode == 4 || mode >= 5) {
          // restore
          if (!weights_ || state_ == -1 || mode == 5) {
               // initialize weights
               delete [] weights_;
               delete alternateWeights_;
               weights_ = new double[numberRows];
               alternateWeights_ = new CoinIndexedVector();
               // enough space so can use it for factorization
               alternateWeights_->reserve(numberRows +
                                          model_->factorization()->maximumPivots());
               if (mode_ != 1 || mode == 5) {
                    // initialize to 1.0 (can we do better?)
                    for (i = 0; i < numberRows; i++) {
                         weights_[i] = 1.0;
                    }
               } else {
                    CoinIndexedVector * temp = new CoinIndexedVector();
                    temp->reserve(numberRows +
                                  model_->factorization()->maximumPivots());
                    double * array = alternateWeights_->denseVector();
                    int * which = alternateWeights_->getIndices();
                    for (i = 0; i < numberRows; i++) {
                         double value = 0.0;
                         array[0] = 1.0;
                         which[0] = i;
                         alternateWeights_->setNumElements(1);
                         alternateWeights_->setPackedMode(true);
                         model_->factorization()->updateColumnTranspose(temp,
                                   alternateWeights_);
                         int number = alternateWeights_->getNumElements();
                         int j;
                         for (j = 0; j < number; j++) {
                              value += array[j] * array[j];
                              array[j] = 0.0;
                         }
                         alternateWeights_->setNumElements(0);
                         weights_[i] = value;
                    }
                    delete temp;
               }
               // create saved weights (not really indexedvector)
               savedWeights_ = new CoinIndexedVector();
               savedWeights_->reserve(numberRows);

               double * array = savedWeights_->denseVector();
               int * which = savedWeights_->getIndices();
               for (i = 0; i < numberRows; i++) {
                    array[i] = weights_[i];
                    which[i] = pivotVariable[i];
               }
          } else if (mode != 6) {
               int * which = alternateWeights_->getIndices();
               CoinIndexedVector * rowArray3 = model_->rowArray(3);
               rowArray3->clear();
               int * back = rowArray3->getIndices();
               // In case something went wrong
               for (i = 0; i < numberRows + numberColumns; i++)
                    back[i] = -1;
               if (mode != 4) {
                    // save
                    CoinMemcpyN(which,	numberRows, savedWeights_->getIndices());
                    CoinMemcpyN(weights_,	numberRows, savedWeights_->denseVector());
               } else {
                    // restore
                    //memcpy(which,savedWeights_->getIndices(),
                    //     numberRows*sizeof(int));
                    //memcpy(weights_,savedWeights_->denseVector(),
                    //     numberRows*sizeof(double));
                    which = savedWeights_->getIndices();
               }
               // restore (a bit slow - but only every re-factorization)
               double * array = savedWeights_->denseVector();
               for (i = 0; i < numberRows; i++) {
                    int iSeq = which[i];
                    back[iSeq] = i;
               }
               for (i = 0; i < numberRows; i++) {
                    int iPivot = pivotVariable[i];
                    iPivot = back[iPivot];
                    if (iPivot >= 0) {
                         weights_[i] = array[iPivot];
                         if (weights_[i] < DEVEX_TRY_NORM)
                              weights_[i] = DEVEX_TRY_NORM; // may need to check more
                    } else {
                         // odd
                         weights_[i] = 1.0;
                    }
               }
          } else {
               // mode 6 - scale back weights as primal errors
               double primalError = model_->largestPrimalError();
               double allowed ;
               if (primalError > 1.0e3)
                    allowed = 10.0;
               else if (primalError > 1.0e2)
                    allowed = 50.0;
               else if (primalError > 1.0e1)
                    allowed = 100.0;
               else
                    allowed = 1000.0;
               double allowedInv = 1.0 / allowed;
               for (i = 0; i < numberRows; i++) {
                    double value = weights_[i];
                    if (value < allowedInv)
                         value = allowedInv;
                    else if (value > allowed)
                         value = allowed;
                    weights_[i] = allowed;
               }
          }
          state_ = 0;
          // set up infeasibilities
          if (!infeasible_) {
               infeasible_ = new CoinIndexedVector();
               infeasible_->reserve(numberRows);
          }
     }
     if (mode >= 2) {
          // Get dubious weights
          //if (!dubiousWeights_)
          //dubiousWeights_=new int[numberRows];
          //model_->factorization()->getWeights(dubiousWeights_);
          infeasible_->clear();
          int iRow;
          const int * pivotVariable = model_->pivotVariable();
          double tolerance = model_->currentPrimalTolerance();
          for (iRow = 0; iRow < numberRows; iRow++) {
               int iPivot = pivotVariable[iRow];
               double value = model_->solution(iPivot);
               double lower = model_->lower(iPivot);
               double upper = model_->upper(iPivot);
               if (value < lower - tolerance) {
                    value -= lower;
                    value *= value;
#ifdef CLP_DUAL_COLUMN_MULTIPLIER
                    if (iPivot < numberColumns)
                         value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
#endif
#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
                    if (lower == upper)
                         value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
#endif
                    // store square in list
                    infeasible_->quickAdd(iRow, value);
               } else if (value > upper + tolerance) {
                    value -= upper;
                    value *= value;
#ifdef CLP_DUAL_COLUMN_MULTIPLIER
                    if (iPivot < numberColumns)
                         value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
#endif
#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
                    if (lower == upper)
                         value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
#endif
                    // store square in list
                    infeasible_->quickAdd(iRow, value);
               }
          }
     }
}
// Gets rid of last update
void
ClpDualRowSteepest::unrollWeights()
{
     double * saved = alternateWeights_->denseVector();
     int number = alternateWeights_->getNumElements();
     int * which = alternateWeights_->getIndices();
     int i;
     if (alternateWeights_->packedMode()) {
          for (i = 0; i < number; i++) {
               int iRow = which[i];
               weights_[iRow] = saved[i];
               saved[i] = 0.0;
          }
     } else {
          for (i = 0; i < number; i++) {
               int iRow = which[i];
               weights_[iRow] = saved[iRow];
               saved[iRow] = 0.0;
          }
     }
     alternateWeights_->setNumElements(0);
}
//-------------------------------------------------------------------
// Clone
//-------------------------------------------------------------------
ClpDualRowPivot * ClpDualRowSteepest::clone(bool CopyData) const
{
     if (CopyData) {
          return new ClpDualRowSteepest(*this);
     } else {
          return new ClpDualRowSteepest();
     }
}
// Gets rid of all arrays
void
ClpDualRowSteepest::clearArrays()
{
     if (persistence_ == normal) {
          delete [] weights_;
          weights_ = NULL;
          delete [] dubiousWeights_;
          dubiousWeights_ = NULL;
          delete infeasible_;
          infeasible_ = NULL;
          delete alternateWeights_;
          alternateWeights_ = NULL;
          delete savedWeights_;
          savedWeights_ = NULL;
     }
     state_ = -1;
}
// Returns true if would not find any row
bool
ClpDualRowSteepest::looksOptimal() const
{
     int iRow;
     const int * pivotVariable = model_->pivotVariable();
     double tolerance = model_->currentPrimalTolerance();
     // we can't really trust infeasibilities if there is primal error
     // this coding has to mimic coding in checkPrimalSolution
     double error = CoinMin(1.0e-2, model_->largestPrimalError());
     // allow tolerance at least slightly bigger than standard
     tolerance = tolerance +  error;
     // But cap
     tolerance = CoinMin(1000.0, tolerance);
     int numberRows = model_->numberRows();
     int numberInfeasible = 0;
     for (iRow = 0; iRow < numberRows; iRow++) {
          int iPivot = pivotVariable[iRow];
          double value = model_->solution(iPivot);
          double lower = model_->lower(iPivot);
          double upper = model_->upper(iPivot);
          if (value < lower - tolerance) {
               numberInfeasible++;
          } else if (value > upper + tolerance) {
               numberInfeasible++;
          }
     }
     return (numberInfeasible == 0);
}
// Called when maximum pivots changes
void
ClpDualRowSteepest::maximumPivotsChanged()
{
     if (alternateWeights_ &&
               alternateWeights_->capacity() != model_->numberRows() +
               model_->factorization()->maximumPivots()) {
          delete alternateWeights_;
          alternateWeights_ = new CoinIndexedVector();
          // enough space so can use it for factorization
          alternateWeights_->reserve(model_->numberRows() +
                                     model_->factorization()->maximumPivots());
     }
}