Changeset 393


Ignore:
Timestamp:
Jul 21, 2004 4:36:05 PM (16 years ago)
Author:
forrest
Message:

Some quadratic stuff

Location:
trunk
Files:
2 added
2 deleted
35 edited

Legend:

Unmodified
Added
Removed
  • trunk/ClpCholeskyWssmp.cpp

    r389 r393  
    488488  //integerParameters_[11]=1;
    489489  integerParameters_[12]=2;
    490   doubleParameters_[10]=max(1.0e-20,largest2);
     490  doubleParameters_[10]=max(1.0e-20,largest);
    491491  if (doubleParameters_[10]>1.0e-3)
    492492    integerParameters_[9]=1;
  • trunk/ClpCholeskyWssmpKKT.cpp

    r389 r393  
    321321    quadratic = quadraticObj->quadraticObjective();
    322322  // matrix
    323   largest=1.0e-100;
    324323  if (!quadratic) {
    325324    for (iColumn=0;iColumn<numberColumns;iColumn++) {
     
    336335          choleskyRow_[numberElements]=row[j]+numberTotal;
    337336          sparseFactor_[numberElements++]=element[j];
     337          largest = max(largest,fabs(element[j]));
    338338        }
    339339      } else {
     
    359359          int jColumn = columnQuadratic[j];
    360360          if (jColumn>iColumn) {
    361             sparseFactor_[numberElements]=quadraticElement[j];
     361            sparseFactor_[numberElements]=-quadraticElement[j];
    362362            choleskyRow_[numberElements++]=jColumn;
    363363          } else if (iColumn==jColumn) {
     
    372372          choleskyRow_[numberElements]=row[j]+numberTotal;
    373373          sparseFactor_[numberElements++]=element[j];
     374          largest = max(largest,fabs(element[j]));
    374375        }
    375376      } else {
     
    416417  integerParameters_[30]=1;
    417418  doubleParameters_[20]=1.0e100;
     419  double largest2= largest*1.0e-20;
     420  largest = min (largest2,1.0e-11);
     421  doubleParameters_[10]=max(1.0e-20,largest);
     422  if (doubleParameters_[10]>1.0e-3)
     423    integerParameters_[9]=1;
     424  else
     425    integerParameters_[9]=0;
    418426#ifndef WSMP
    419427  // Set up LDL cutoff
    420428  integerParameters_[34]=numberTotal;
    421   doubleParameters_[10]=min(largest*1.0e-19,1.0e-14);
    422429  doubleParameters_[20]=1.0e-15;
    423430  doubleParameters_[34]=1.0e-12;
  • trunk/ClpInterior.cpp

    r389 r393  
    8888  deltaY_(NULL),
    8989  deltaZ_(NULL),
    90   deltaT_(NULL),
     90  deltaW_(NULL),
    9191  deltaSU_(NULL),
    9292  deltaSL_(NULL),
     
    9797  rhsL_(NULL),
    9898  rhsZ_(NULL),
    99   rhsT_(NULL),
     99  rhsW_(NULL),
    100100  rhsC_(NULL),
    101101  zVec_(NULL),
    102   tVec_(NULL),
     102  wVec_(NULL),
    103103  cholesky_(NULL),
    104104  numberComplementarityPairs_(0),
     
    177177    deltaY_(NULL),
    178178    deltaZ_(NULL),
    179     deltaT_(NULL),
     179    deltaW_(NULL),
    180180    deltaSU_(NULL),
    181181    deltaSL_(NULL),
     
    186186    rhsL_(NULL),
    187187    rhsZ_(NULL),
    188     rhsT_(NULL),
     188    rhsW_(NULL),
    189189    rhsC_(NULL),
    190190    zVec_(NULL),
    191     tVec_(NULL),
     191    wVec_(NULL),
    192192    cholesky_(NULL),
    193193    numberComplementarityPairs_(0),
     
    252252  deltaY_(NULL),
    253253  deltaZ_(NULL),
    254   deltaT_(NULL),
     254  deltaW_(NULL),
    255255  deltaSU_(NULL),
    256256  deltaSL_(NULL),
     
    261261  rhsL_(NULL),
    262262  rhsZ_(NULL),
    263   rhsT_(NULL),
     263  rhsW_(NULL),
    264264  rhsC_(NULL),
    265265  zVec_(NULL),
    266   tVec_(NULL),
     266  wVec_(NULL),
    267267  cholesky_(NULL)
    268268{
     
    330330  deltaY_(NULL),
    331331  deltaZ_(NULL),
    332   deltaT_(NULL),
     332  deltaW_(NULL),
    333333  deltaSU_(NULL),
    334334  deltaSL_(NULL),
     
    339339  rhsL_(NULL),
    340340  rhsZ_(NULL),
    341   rhsT_(NULL),
     341  rhsW_(NULL),
    342342  rhsC_(NULL),
    343343  zVec_(NULL),
    344   tVec_(NULL),
     344  wVec_(NULL),
    345345  cholesky_(NULL),
    346346  numberComplementarityPairs_(0),
     
    427427  deltaX_ = ClpCopyOfArray(rhs.deltaX_,numberRows_+numberColumns_);
    428428  deltaZ_ = ClpCopyOfArray(rhs.deltaZ_,numberRows_+numberColumns_);
    429   deltaT_ = ClpCopyOfArray(rhs.deltaT_,numberRows_+numberColumns_);
     429  deltaW_ = ClpCopyOfArray(rhs.deltaW_,numberRows_+numberColumns_);
    430430  deltaSU_ = ClpCopyOfArray(rhs.deltaSU_,numberRows_+numberColumns_);
    431431  deltaSL_ = ClpCopyOfArray(rhs.deltaSL_,numberRows_+numberColumns_);
     
    436436  rhsL_ = ClpCopyOfArray(rhs.rhsL_,numberRows_+numberColumns_);
    437437  rhsZ_ = ClpCopyOfArray(rhs.rhsZ_,numberRows_+numberColumns_);
    438   rhsT_ = ClpCopyOfArray(rhs.rhsT_,numberRows_+numberColumns_);
     438  rhsW_ = ClpCopyOfArray(rhs.rhsW_,numberRows_+numberColumns_);
    439439  rhsC_ = ClpCopyOfArray(rhs.rhsC_,numberRows_+numberColumns_);
    440440  solution_ = ClpCopyOfArray(rhs.solution_,numberRows_+numberColumns_);
    441441  workArray_ = ClpCopyOfArray(rhs.workArray_,numberRows_+numberColumns_);
    442442  zVec_ = ClpCopyOfArray(rhs.zVec_,numberRows_+numberColumns_);
    443   tVec_ = ClpCopyOfArray(rhs.tVec_,numberRows_+numberColumns_);
     443  wVec_ = ClpCopyOfArray(rhs.wVec_,numberRows_+numberColumns_);
    444444  cholesky_ = rhs.cholesky_->clone();
    445445  numberComplementarityPairs_ = rhs.numberComplementarityPairs_;
     
    494494  delete [] deltaZ_;
    495495  deltaZ_ = NULL;
    496   delete [] deltaT_;
    497   deltaT_ = NULL;
     496  delete [] deltaW_;
     497  deltaW_ = NULL;
    498498  delete [] deltaSU_;
    499499  deltaSU_ = NULL;
     
    512512  delete [] rhsZ_;
    513513  rhsZ_ = NULL;
    514   delete [] rhsT_;
    515   rhsT_ = NULL;
     514  delete [] rhsW_;
     515  rhsW_ = NULL;
    516516  delete [] rhsC_;
    517517  rhsC_ = NULL;
     
    522522  delete [] zVec_;
    523523  zVec_ = NULL;
    524   delete [] tVec_;
    525   tVec_ = NULL;
     524  delete [] wVec_;
     525  wVec_ = NULL;
    526526  delete cholesky_;
    527527}
     
    632632  deltaZ_ = new double [nTotal];
    633633  CoinZeroN(deltaZ_,nTotal);
    634   assert (!deltaT_);
    635   deltaT_ = new double [nTotal];
    636   CoinZeroN(deltaT_,nTotal);
     634  assert (!deltaW_);
     635  deltaW_ = new double [nTotal];
     636  CoinZeroN(deltaW_,nTotal);
    637637  assert (!deltaSU_);
    638638  deltaSU_ = new double [nTotal];
     
    662662  rhsZ_ = new double [nTotal];
    663663  CoinZeroN(rhsZ_,nTotal);
    664   assert (!rhsT_);
    665   rhsT_ = new double [nTotal];
    666   CoinZeroN(rhsT_,nTotal);
     664  assert (!rhsW_);
     665  rhsW_ = new double [nTotal];
     666  CoinZeroN(rhsW_,nTotal);
    667667  assert (!rhsC_);
    668668  rhsC_ = new double [nTotal];
     
    674674  zVec_ = new double [nTotal];
    675675  CoinZeroN(zVec_,nTotal);
    676   assert (!tVec_);
    677   tVec_ = new double [nTotal];
    678   CoinZeroN(tVec_,nTotal);
     676  assert (!wVec_);
     677  wVec_ = new double [nTotal];
     678  CoinZeroN(wVec_,nTotal);
    679679  assert (!dj_);
    680680  dj_ = new double [nTotal];
     
    749749  delete [] zVec_;
    750750  zVec_ = NULL;
    751   delete [] tVec_;
    752   tVec_ = NULL;
     751  delete [] wVec_;
     752  wVec_ = NULL;
    753753  delete [] dj_;
    754754  dj_ = NULL;
  • trunk/ClpLinearObjective.cpp

    r225 r393  
    33
    44#include "CoinPragma.hpp"
    5 #include "ClpModel.hpp"
     5#include "CoinIndexedVector.hpp"
     6#include "ClpFactorization.hpp"
     7#include "ClpSimplex.hpp"
    68#include "ClpLinearObjective.hpp"
    79
     
    110112// Returns gradient
    111113double * 
    112 ClpLinearObjective::gradient(const double * solution, double & offset,bool refresh)
    113 {
     114ClpLinearObjective::gradient(const ClpSimplex * model,
     115                             const double * solution, double & offset,bool refresh,
     116                             int includeLinear)
     117{
     118  // not sure what to do about scaling
     119  //assert (!model);
     120  assert (includeLinear==2); //otherwise need to return all zeros
    114121  offset=0.0;
    115122  return objective_;
    116123}
    117124 
     125/* Returns reduced gradient.Returns an offset (to be added to current one).
     126 */
     127double
     128ClpLinearObjective::reducedGradient(ClpSimplex * model, double * region,
     129                                 bool useFeasibleCosts)
     130{
     131  int numberRows = model->numberRows();
     132  //work space
     133  CoinIndexedVector  * workSpace = model->rowArray(0);
     134 
     135  CoinIndexedVector arrayVector;
     136  arrayVector.reserve(numberRows+1);
     137 
     138  int iRow;
     139#ifdef CLP_DEBUG
     140  workSpace->checkClear();
     141#endif
     142  double * array = arrayVector.denseVector();
     143  int * index = arrayVector.getIndices();
     144  int number=0;
     145  const double * cost = model->costRegion();
     146  assert (!useFeasibleCosts);
     147  const int * pivotVariable = model->pivotVariable();
     148  for (iRow=0;iRow<numberRows;iRow++) {
     149    int iPivot=pivotVariable[iRow];
     150    double value = cost[iPivot];
     151    if (value) {
     152      array[iRow]=value;
     153      index[number++]=iRow;
     154    }
     155  }
     156  arrayVector.setNumElements(number);
     157
     158  int numberColumns=model->numberColumns();
     159 
     160  // Btran basic costs
     161  double * work = workSpace->denseVector();
     162  model->factorization()->updateColumnTranspose(workSpace,&arrayVector);
     163  ClpFillN(work,numberRows,0.0);
     164  // now look at dual solution
     165  double * rowReducedCost = region+numberColumns;
     166  double * dual = rowReducedCost;
     167  double * rowCost = model->costRegion(0);
     168  for (iRow=0;iRow<numberRows;iRow++) {
     169    dual[iRow]=array[iRow];
     170  }
     171  double * dj = region;
     172  ClpDisjointCopyN(model->costRegion(1),numberColumns,dj);
     173  model->transposeTimes(-1.0,dual,dj);
     174  for (iRow=0;iRow<numberRows;iRow++) {
     175    // slack
     176    double value = dual[iRow];
     177    value += rowCost[iRow];
     178    rowReducedCost[iRow]=value;
     179  }
     180  return 0.0;
     181}
     182/* Returns step length which gives minimum of objective for
     183   solution + theta * change vector up to maximum theta.
     184   
     185   arrays are numberColumns+numberRows
     186*/
     187double
     188ClpLinearObjective::stepLength(ClpSimplex * model,
     189                               const double * solution,
     190                               const double * change,
     191                               double maximumTheta,
     192                               double & currentObj,
     193                               double & predictedObj,
     194                               double & thetaObj)
     195{
     196  const double * cost = model->costRegion();
     197  double delta=0.0;
     198  int numberRows = model->numberRows();
     199  int numberColumns = model->numberColumns();
     200  currentObj=0.0;
     201  thetaObj=0.0;
     202  for (int iColumn=0;iColumn<numberColumns+numberRows;iColumn++) {
     203    delta += cost[iColumn]*change[iColumn];
     204    currentObj += cost[iColumn]*solution[iColumn];
     205  }
     206  thetaObj =currentObj + delta*maximumTheta;
     207  predictedObj =currentObj + delta*maximumTheta;
     208  if (delta<0.0) {
     209    return maximumTheta;
     210  } else {
     211    printf("odd linear direction %g\n",delta);
     212    return 0.0;
     213  }
     214}
    118215//-------------------------------------------------------------------
    119216// Clone
     
    180277  }
    181278}
    182 
     279// Scale objective
     280void
     281ClpLinearObjective::reallyScale(const double * columnScale)
     282{
     283  for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
     284    objective_[iColumn] *= columnScale[iColumn];
     285  }
     286}
     287
  • trunk/ClpMatrixBase.cpp

    r336 r393  
    405405{
    406406}
    407 
    408 
     407// Really scale matrix
     408void
     409ClpMatrixBase::reallyScale(const double * rowScale, const double * columnScale)
     410{
     411  std::cerr<<"reallyScale not supported - ClpMatrixBase"<<std::endl;
     412  abort();
     413}
     414
     415
  • trunk/ClpMessage.cpp

    r371 r393  
    9595  {CLP_BARRIER_STEP,58,2,"Steps - primal %g ,dual %g , mu %g"},
    9696  {CLP_RIM_SCALE,59,1,"Automatic rim scaling gives objective scale of %g and rhs/bounds scale of %g"},
     97  {CLP_SLP_ITER,58,1,"Pass %d objective %g - drop %g, largest delta %g"},
    9798  {CLP_DUMMY_END,999999,0,""}
    9899};
  • trunk/ClpModel.cpp

    r389 r393  
    286286//#############################################################################
    287287// Copy constructor.
    288 ClpModel::ClpModel(const ClpModel &rhs) :
     288ClpModel::ClpModel(const ClpModel &rhs, int scalingMode) :
    289289  optimizationDirection_(rhs.optimizationDirection_),
    290290  numberRows_(rhs.numberRows_),
     
    292292{
    293293  gutsOfCopy(rhs);
     294  if (scalingMode>=0&&matrix_&&
     295      matrix_->allElementsInRange(this,smallElement_,1.0e20)) {
     296    // really do scaling
     297    scalingFlag_=scalingMode;
     298    delete [] rowScale_;
     299    rowScale_ = NULL;
     300    delete [] columnScale_;
     301    columnScale_ = NULL;
     302    if (!matrix_->scale(this)) {
     303      // scaling worked - now apply
     304      gutsOfScaling();
     305      // pretend not scaled
     306      scalingFlag_ = -scalingFlag_;
     307    } else {
     308      // not scaled
     309      scalingFlag_=0;
     310    }
     311  }
    294312  CoinSeedRandom(1234567);
    295313}
     
    417435    rowCopy_ = NULL;
    418436    ray_ = rhs.ray_;
    419     rowScale_ = rhs.rowScale_;
    420     columnScale_ = rhs.columnScale_;
     437    //rowScale_ = rhs.rowScale_;
     438    //columnScale_ = rhs.columnScale_;
    421439    lengthNames_ = 0;
    422440    rowNames_ = std::vector<std::string> ();
     
    467485  otherModel.ray_ = ray_;
    468486  ray_ = NULL;
    469   rowScale_=NULL;
    470   columnScale_=NULL;
     487  //rowScale_=NULL;
     488  //columnScale_=NULL;
    471489  // do status
    472490  if (otherModel.status_!=status_) {
     
    13491367  assert ((dynamic_cast< ClpLinearObjective*>(objective_)));
    13501368  double offset;
    1351   ClpObjective * obj = new ClpQuadraticObjective(objective_->gradient(NULL,offset,false),numberColumns,
    1352                                              start,column,element);
     1369  ClpObjective * obj = new ClpQuadraticObjective(objective_->gradient(NULL,NULL,offset,false),
     1370                                                 numberColumns,
     1371                                                 start,column,element);
    13531372  delete objective_;
    13541373  objective_ = obj;
     
    13621381  double offset;
    13631382  ClpQuadraticObjective * obj =
    1364     new ClpQuadraticObjective(objective_->gradient(NULL,offset,false),numberColumns_,
    1365                                                  NULL,NULL,NULL);
     1383    new ClpQuadraticObjective(objective_->gradient(NULL,NULL,offset,false),
     1384                              numberColumns_,
     1385                              NULL,NULL,NULL);
    13661386  delete objective_;
    13671387  objective_ = obj;
     
    18301850  writer.setObjectiveOffset(objectiveOffset());
    18311851  delete [] objective;
    1832   return writer.writeMps(filename, 0 /* do not gzip it*/, formatType, numberAcross);
     1852  // allow for quadratic objective
     1853#ifndef NO_RTTI
     1854  ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
     1855#else
     1856  ClpQuadraticObjective * quadraticObj = NULL;
     1857  if (objective_->type()==2)
     1858    quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
     1859#endif
     1860  CoinPackedMatrix * quadratic=NULL;
     1861  if (quadraticObj)
     1862    quadratic = quadraticObj->quadraticObjective();
     1863  return writer.writeMps(filename, 0 /* do not gzip it*/, formatType, numberAcross,
     1864                         quadratic);
    18331865  if (rowNames) {
    18341866    for (int iRow=0;iRow<numberRows_;iRow++) {
     
    18811913    matrix_->transposeTimes(scalar,x,y);
    18821914}
     1915// Does much of scaling
     1916void
     1917ClpModel::gutsOfScaling()
     1918{
     1919  int i;
     1920  if (rowObjective_) {
     1921    for (i=0;i<numberRows_;i++)
     1922      rowObjective_[i] /= rowScale_[i];
     1923  }
     1924  for (i=0;i<numberRows_;i++) {
     1925    double multiplier = rowScale_[i];
     1926    double inverseMultiplier = 1.0/multiplier;
     1927    rowActivity_[i] *= multiplier;
     1928    dual_[i] *= inverseMultiplier;
     1929    if (rowLower_[i]>-1.0e30)
     1930      rowLower_[i] *= multiplier;
     1931    else
     1932      rowLower_[i] = -COIN_DBL_MAX;
     1933    if (rowUpper_[i]<1.0e30)
     1934      rowUpper_[i] *= multiplier;
     1935    else
     1936      rowUpper_[i] = COIN_DBL_MAX;
     1937  }
     1938  for (i=0;i<numberColumns_;i++) {
     1939    double multiplier = 1.0/columnScale_[i];
     1940    columnActivity_[i] *= multiplier;
     1941    reducedCost_[i] *= columnScale_[i];
     1942    if (columnLower_[i]>-1.0e30)
     1943      columnLower_[i] *= multiplier;
     1944    else
     1945      columnLower_[i] = -COIN_DBL_MAX;
     1946    if (columnUpper_[i]<1.0e30)
     1947      columnUpper_[i] *= multiplier;
     1948    else
     1949      columnUpper_[i] = COIN_DBL_MAX;
     1950   
     1951  }
     1952  //now replace matrix
     1953  //and objective
     1954  matrix_->reallyScale(rowScale_,columnScale_);
     1955  objective_->reallyScale(columnScale_);
     1956}
     1957/* If we constructed a "really" scaled model then this reverses the operation.
     1958      Quantities may not be exactly as they were before due to rounding errors */
     1959void
     1960ClpModel::unscale()
     1961{
     1962  if (rowScale_) {
     1963    int i;
     1964    // reverse scaling
     1965    for (i=0;i<numberRows_;i++)
     1966      rowScale_[i] = 1.0/rowScale_[i];
     1967    for (i=0;i<numberColumns_;i++)
     1968      columnScale_[i] = 1.0/columnScale_[i];
     1969    gutsOfScaling();
     1970  }
     1971 
     1972  scalingFlag_=0;
     1973  delete [] rowScale_;
     1974  rowScale_ = NULL;
     1975  delete [] columnScale_;
     1976  columnScale_ = NULL;
     1977}
  • trunk/ClpNonLinearCost.cpp

    r336 r393  
    455455          if (value<lowerValue-primalTolerance) {
    456456            value = lowerValue-value;
     457#ifndef NDEBUG
     458            if(value>1.0e15)
     459              printf("nonlincostb %d %g %g %g\n",
     460                     iSequence,lowerValue,solution[iSequence],lower_[iRange+2]);
     461#endif
    457462            sumInfeasibilities_ += value;
    458463            largestInfeasibility_ = max(largestInfeasibility_,value);
     
    467472          if (value>upperValue+primalTolerance) {
    468473            value = value-upperValue;
     474#ifndef NDEBUG
     475            if(value>1.0e15)
     476              printf("nonlincostu %d %g %g %g\n",
     477                     iSequence,lower_[iRange-1],solution[iSequence],upperValue);
     478#endif
    469479            sumInfeasibilities_ += value;
    470480            largestInfeasibility_ = max(largestInfeasibility_,value);
     
    985995  double value;
    986996  model_->getDblParam(ClpObjOffset,value);
    987   return feasibleCost_*model_->optimizationDirection()/
     997  return (feasibleCost_+model_->objectiveAsObject()->nonlinearOffset())*model_->optimizationDirection()/
    988998    (model_->objectiveScale()*model_->rhsScale())-value;
    989999}
  • trunk/ClpObjective.cpp

    r240 r393  
    1414//-------------------------------------------------------------------
    1515ClpObjective::ClpObjective () :
    16   type_(-1)
     16  offset_(0.0),
     17  type_(-1),
     18  activated_(1)
    1719{
    1820
     
    2325//-------------------------------------------------------------------
    2426ClpObjective::ClpObjective (const ClpObjective & source) :
    25   type_(source.type_)
     27  offset_(source.offset_),
     28  type_(source.type_),
     29  activated_(source.activated_)
    2630
    2731
     
    4347{
    4448  if (this != &rhs) {
     49    offset_ = rhs.offset_;
    4550    type_ = rhs.type_;
     51    activated_= rhs.activated_;
    4652  }
    4753  return *this;
     
    5864  return NULL;
    5965}
     66/* Given a zeroed array sets nonlinear columns to 1.
     67   Returns number of nonlinear columns
     68*/
     69int
     70ClpObjective::markNonlinear(char * which)
     71{
     72  return 0;
     73}
    6074
  • trunk/ClpPackedMatrix.cpp

    r372 r393  
    22572257  return copy;
    22582258}
     2259// Really scale matrix
     2260void
     2261ClpPackedMatrix::reallyScale(const double * rowScale, const double * columnScale)
     2262{
     2263  int numberColumns = matrix_->getNumCols();
     2264  const int * row = matrix_->getIndices();
     2265  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
     2266  const int * length = matrix_->getVectorLengths();
     2267  double * element = matrix_->getMutableElements();
     2268  // scale
     2269  for (int iColumn=0;iColumn<numberColumns;iColumn++) {
     2270    CoinBigIndex j;
     2271    double scale = columnScale[iColumn];
     2272    for (j=columnStart[iColumn];j<columnStart[iColumn]+length[iColumn];j++) {
     2273      int iRow = row[j];
     2274      element[j] *= scale*rowScale[iRow];
     2275    }
     2276  }
     2277}
  • trunk/ClpPlusMinusOneMatrix.cpp

    r360 r393  
    869869  }
    870870  columnArray->setNumElements(numberNonZero);
     871  if (packed)
     872    columnArray->setPacked();
    871873  y->setNumElements(0);
    872874}
  • trunk/ClpPredictorCorrector.cpp

    r389 r393  
    6161  int modeSwitch=0;
    6262  //if (quadraticObj)
    63   //modeSwitch |= 1; // switch off centring for now
     63    //modeSwitch |= 1; // switch off centring for now
    6464  //if (quadraticObj)
    6565  //modeSwitch |=4;
     
    152152  double * saveY = new double[numberRows_];
    153153  double * saveZ = new double[numberTotal];
    154   double * saveT = new double[numberTotal];
     154  double * saveW = new double[numberTotal];
    155155  double * saveSL = new double[numberTotal];
    156156  double * saveSU = new double[numberTotal];
     
    170170      printf(" col  dsol  ddiag  dwvec  dzvec dbdslu dbdsll\n");
    171171      for (i=0;i<numberColumns_+numberRows_;i++) {
    172         printf(" %d %g %g %g %g %g %g\n",i,solution_[i],diagonal_[i],tVec_[i],
     172        printf(" %d %g %g %g %g %g %g\n",i,solution_[i],diagonal_[i],wVec_[i],
    173173               zVec_[i],upperSlack_[i],lowerSlack_[i]);
    174174      }
     
    281281        <<CoinMessageEol;
    282282    }
    283     //if (complementarityGap_>=0.98*lastComplementarityGap) {
     283    int numberBack = quadraticObj ? 10 : 5;
    284284    //tryJustPredictor=true;
    285285    //printf("trying just predictor\n");
     
    294294          <<CoinMessageEol;
    295295        break;
    296       } else if (numberIterations_-lastGoodIteration>=5&&
     296      } else if (numberIterations_-lastGoodIteration>=numberBack&&
    297297                 complementarityGap_<1.0e-6) {
    298298        break; // not doing very well - give up
     
    301301      lastGoodIteration=numberIterations_;
    302302      lastComplementarityGap=complementarityGap_;
    303     } else if (numberIterations_-lastGoodIteration>=5&&
     303    } else if (numberIterations_-lastGoodIteration>=numberBack&&
    304304               complementarityGap_<1.0e-3) {
    305305      handler_->message(CLP_BARRIER_COMPLEMENTARITY,messages_)
     
    421421    //set up for affine direction
    422422    setupForSolve(phase);
    423     directionAccuracy=findDirectionVector(phase);
    424     if (directionAccuracy>worstDirectionAccuracy_) {
    425       worstDirectionAccuracy_=directionAccuracy;
    426     }
    427     if (saveIteration>0&&directionAccuracy>1.0) {
    428       handler_->message(CLP_BARRIER_EXIT2,messages_)
    429         <<saveIteration
    430         <<CoinMessageEol;
    431       break;
    432     }
    433     findStepLength(phase);
    434     nextGap=complementarityGap(nextNumber,nextNumberItems,1);
    435     debugMove(0,actualPrimalStep_,actualDualStep_);
    436     //debugMove(0,1.0e-7,1.0e-7);
     423    if ((modeSwitch&2)==0) {
     424      directionAccuracy=findDirectionVector(phase);
     425      if (directionAccuracy>worstDirectionAccuracy_) {
     426        worstDirectionAccuracy_=directionAccuracy;
     427      }
     428      if (saveIteration>0&&directionAccuracy>1.0) {
     429        handler_->message(CLP_BARRIER_EXIT2,messages_)
     430          <<saveIteration
     431          <<CoinMessageEol;
     432        break;
     433      }
     434      findStepLength(phase);
     435      nextGap=complementarityGap(nextNumber,nextNumberItems,1);
     436      debugMove(0,actualPrimalStep_,actualDualStep_);
     437      debugMove(0,1.0e-2,1.0e-2);
     438    }
    437439    double affineGap=nextGap;
    438440    int bestPhase=0;
     
    440442    // ?
    441443    bestNextGap=max(nextGap,0.8*complementarityGap_);
    442     //bestNextGap=max(nextGap,0.99*complementarityGap_);
     444    if (quadraticObj)
     445      bestNextGap=max(nextGap,0.99*complementarityGap_);
    443446    if (complementarityGap_>1.0e-4*numberComplementarityPairs_) {
    444447      //std::cout <<"predicted duality gap "<<nextGap<<std::endl;
     
    520523    }
    521524    // If bad gap - try standard primal dual
    522     if (nextGap>complementarityGap_*1.0001)
     525    if (nextGap>complementarityGap_*1.001)
    523526      goodMove=false;
    524527    if ((modeSwitch&2)!=0)
     
    528531      memcpy(saveY,deltaY_,numberRows_*sizeof(double));
    529532      memcpy(saveZ,deltaZ_,numberTotal*sizeof(double));
    530       memcpy(saveT,deltaT_,numberTotal*sizeof(double));
     533      memcpy(saveW,deltaW_,numberTotal*sizeof(double));
    531534      memcpy(saveSL,deltaSL_,numberTotal*sizeof(double));
    532535      memcpy(saveSU,deltaSU_,numberTotal*sizeof(double));
     
    561564      if (goodMove) {
    562565        phase=1;
    563         findStepLength(phase);
     566        double norm = findStepLength(phase);
    564567        nextGap = complementarityGap(nextNumber,nextNumberItems,1);
    565568        debugMove(1,actualPrimalStep_,actualDualStep_);
    566569        //debugMove(1,1.0e-7,1.0e-7);
    567570        goodMove=checkGoodMove(true,bestNextGap,allowIncreasingGap);
     571        if (norm<0)
     572          goodMove=false;
    568573        if (!goodMove) {
    569574#ifdef SOME_DEBUG
     
    589594          deltaX_[i] = lambda*deltaX_[i]+(1.0-lambda)*saveX[i];
    590595          deltaZ_[i] = lambda*deltaZ_[i]+(1.0-lambda)*saveZ[i];
    591           deltaT_[i] = lambda*deltaT_[i]+(1.0-lambda)*saveT[i];
     596          deltaW_[i] = lambda*deltaW_[i]+(1.0-lambda)*saveW[i];
    592597          deltaSL_[i] = lambda*deltaSL_[i]+(1.0-lambda)*saveSL[i];
    593598          deltaSU_[i] = lambda*deltaSU_[i]+(1.0-lambda)*saveSU[i];
     
    596601        //memcpy(deltaY_,saveY,numberRows_*sizeof(double));
    597602        //memcpy(deltaZ_,saveZ,numberTotal*sizeof(double));
    598         //memcpy(deltaT_,saveT,numberTotal*sizeof(double));
     603        //memcpy(deltaW_,saveW,numberTotal*sizeof(double));
    599604        //memcpy(deltaSL_,saveSL,numberTotal*sizeof(double));
    600605        //memcpy(deltaSU_,saveSU,numberTotal*sizeof(double));
     
    640645      //mu_=min(smallestPrimalDualMu*0.95,mu_);
    641646      smallestPrimalDualMu = mu_;
     647      // Try simpler
     648      floatNumber = numberComplementarityItems_;
     649      mu_=0.5*complementarityGap_/floatNumber;
     650      //if ((modeSwitch&2)==0) {
     651      //if ((numberIterations_&1)==0)
     652      //  mu_ *= 0.5;
     653      //} else {
     654      //mu_ *= 0.8;
     655      //}
    642656      //set up for next step
    643657      setupForSolve(2);
    644658      findDirectionVector(2);
    645       findStepLength(2);
     659      double norm=findStepLength(2);
    646660      // just for debug
    647       bestNextGap = complementarityGap_;
     661      bestNextGap = complementarityGap_*1.0005;
     662      //bestNextGap=COIN_DBL_MAX;
    648663      nextGap=complementarityGap(nextNumber,nextNumberItems,2);
    649664      debugMove(2,actualPrimalStep_,actualDualStep_);
    650665      //debugMove(2,1.0e-7,1.0e-7);
    651666      checkGoodMove(false, bestNextGap,allowIncreasingGap);
    652       if (nextGap>0.9*complementarityGap_&&bestPhase==0&&affineGap<nextGap
    653           &&numberIterations_>80) {
     667      if ((nextGap>0.9*complementarityGap_&&bestPhase==0&&affineGap<nextGap
     668          &&(numberIterations_>80||(numberIterations_>20&&quadraticObj)))||norm<0.0) {
    654669        // Back to affine
    655670        phase=0;
     
    666681    if (!goodMove)
    667682      mu_=nextGap / ((double) 1.1*nextNumber);
    668     //goodMove=true; // don't do if was bad move
     683    //if (quadraticObj)
     684    //goodMove=true;
    669685    //goodMove=false; //TEMP
    670686    // Do centering steps
     
    684700      memcpy(saveY,deltaY_,numberRows_*sizeof(double));
    685701      memcpy(saveZ,deltaZ_,numberTotal*sizeof(double));
    686       memcpy(saveT,deltaT_,numberTotal*sizeof(double));
     702      memcpy(saveW,deltaW_,numberTotal*sizeof(double));
    687703      double savePrimalStep = actualPrimalStep_;
    688704      double saveDualStep = actualDualStep_;
     
    717733        memcpy(deltaY_,saveY,numberRows_*sizeof(double));
    718734        memcpy(deltaZ_,saveZ,numberTotal*sizeof(double));
    719         memcpy(deltaT_,saveT,numberTotal*sizeof(double));
     735        memcpy(deltaW_,saveW,numberTotal*sizeof(double));
    720736      } else {
    721737#ifdef SOME_DEBUG
     
    748764  delete [] saveY;
    749765  delete [] saveZ;
    750   delete [] saveT;
     766  delete [] saveW;
    751767  delete [] saveSL;
    752768  delete [] saveSU;
     
    810826  int chosenPrimalSequence=-1;
    811827  int chosenDualSequence=-1;
    812   double * zVec = zVec_;
    813   double * tVec = tVec_;
    814   double * lowerSlack = lowerSlack_;
    815   double * upperSlack = upperSlack_;
    816   //direction vector in deltaX
    817   double * deltaX = deltaX_;
     828  bool lowPrimal=false;
     829  bool lowDual=false;
    818830  // If done many iterations then allow to hit boundary
    819831  double hitTolerance;
     
    824836    hitTolerance = max(1.0e3,1.0e-3*objectiveNorm_);
    825837  int iColumn;
     838  //printf("dual value %g\n",dual_[0]);
     839  //printf("     X     dX      lS     dlS     uS     dUs    dj    Z dZ     t   dT\n");
    826840  for (iColumn=0;iColumn<numberTotal;iColumn++) {
    827841    if (!flagged(iColumn)) {
    828       double directionElement=deltaX[iColumn];
     842      double directionElement=deltaX_[iColumn];
    829843      if (directionNorm<fabs(directionElement)) {
    830844        directionNorm=fabs(directionElement);
    831       }
     845      }
     846      if (0)
     847      printf("%d %g %g %g %g %g %g %g %g %g %g %g\n",
     848             iColumn,solution_[iColumn],deltaX_[iColumn],
     849             lowerSlack_[iColumn],deltaSL_[iColumn],
     850             upperSlack_[iColumn],deltaSU_[iColumn],
     851             dj_[iColumn],
     852             zVec_[iColumn],deltaZ_[iColumn],
     853             wVec_[iColumn],deltaW_[iColumn]);
    832854      if (lowerBound(iColumn)) {
    833855        double delta = - deltaSL_[iColumn];
    834856        double z1 = deltaZ_[iColumn];
    835         double newZ = zVec[iColumn]+z1;
    836         if (zVec[iColumn]>tolerance) {
    837           if (zVec[iColumn]<-z1*maximumDualStep) {
    838             maximumDualStep=-zVec[iColumn]/z1;
     857        double newZ = zVec_[iColumn]+z1;
     858        if (zVec_[iColumn]>tolerance) {
     859          if (zVec_[iColumn]<-z1*maximumDualStep) {
     860            maximumDualStep=-zVec_[iColumn]/z1;
    839861            chosenDualSequence=iColumn;
     862            lowDual=true;
    840863          }
    841864        }
    842         if (lowerSlack[iColumn]<maximumPrimalStep*delta) {
    843           double newStep=lowerSlack[iColumn]/delta;
     865        if (lowerSlack_[iColumn]<maximumPrimalStep*delta) {
     866          double newStep=lowerSlack_[iColumn]/delta;
    844867          if (newStep>0.2||newZ<hitTolerance||delta>1.0e3||delta<=1.0e-6||dj_[iColumn]<hitTolerance) {
    845868            maximumPrimalStep = newStep;
    846869            chosenPrimalSequence=iColumn;
     870            lowPrimal=true;
    847871          } else {
    848872            //printf("small %d delta %g newZ %g step %g\n",iColumn,delta,newZ,newStep);
     
    852876      if (upperBound(iColumn)) {
    853877        double delta = - deltaSU_[iColumn];;
    854         double t1 = deltaT_[iColumn];
    855         double newT = tVec[iColumn]+t1;
    856         if (tVec[iColumn]>tolerance) {
    857           if (tVec[iColumn]<-t1*maximumDualStep) {
    858             maximumDualStep=-tVec[iColumn]/t1;
     878        double w1 = deltaW_[iColumn];
     879        double newT = wVec_[iColumn]+w1;
     880        if (wVec_[iColumn]>tolerance) {
     881          if (wVec_[iColumn]<-w1*maximumDualStep) {
     882            maximumDualStep=-wVec_[iColumn]/w1;
    859883            chosenDualSequence=iColumn;
     884            lowDual=false;
    860885          }
    861886        }
    862         if (upperSlack[iColumn]<maximumPrimalStep*delta) {
    863           double newStep=upperSlack[iColumn]/delta;
     887        if (upperSlack_[iColumn]<maximumPrimalStep*delta) {
     888          double newStep=upperSlack_[iColumn]/delta;
    864889          if (newStep>0.2||newT<hitTolerance||delta>1.0e3||delta<=1.0e-6||dj_[iColumn]>-hitTolerance) {
    865890            maximumPrimalStep = newStep;
    866891            chosenPrimalSequence=iColumn;
     892            lowPrimal=false;
    867893          } else {
    868894            //printf("small %d delta %g newT %g step %g\n",iColumn,delta,newT,newStep);
     
    875901  printf("new step - phase %d, norm %.18g, dual step %.18g, primal step %.18g\n",
    876902         phase,directionNorm,maximumDualStep,maximumPrimalStep);
     903  if (lowDual)
     904    printf("ld %d %g %g => %g (dj %g,sol %g) ",
     905           chosenDualSequence,zVec_[chosenDualSequence],
     906           deltaZ_[chosenDualSequence],zVec_[chosenDualSequence]+
     907           maximumDualStep*deltaZ_[chosenDualSequence],dj_[chosenDualSequence],
     908           solution_[chosenDualSequence]);
     909  else
     910    printf("ud %d %g %g => %g (dj %g,sol %g) ",
     911           chosenDualSequence,wVec_[chosenDualSequence],
     912           deltaW_[chosenDualSequence],wVec_[chosenDualSequence]+
     913           maximumDualStep*deltaW_[chosenDualSequence],dj_[chosenDualSequence],
     914           solution_[chosenDualSequence]);
     915  if (lowPrimal)
     916    printf("lp %d %g %g => %g (dj %g,sol %g)\n",
     917           chosenPrimalSequence,lowerSlack_[chosenPrimalSequence],
     918           deltaSL_[chosenPrimalSequence],lowerSlack_[chosenPrimalSequence]+
     919           maximumPrimalStep*deltaSL_[chosenPrimalSequence],
     920           dj_[chosenPrimalSequence],solution_[chosenPrimalSequence]);
     921  else
     922    printf("up %d %g %g => %g (dj %g,sol %g)\n",
     923           chosenPrimalSequence,upperSlack_[chosenPrimalSequence],
     924           deltaSU_[chosenPrimalSequence],upperSlack_[chosenPrimalSequence]+
     925           maximumPrimalStep*deltaSU_[chosenPrimalSequence],
     926           dj_[chosenPrimalSequence],solution_[chosenPrimalSequence]);
    877927#endif
    878928  actualPrimalStep_=stepLength_*maximumPrimalStep;
     
    893943#endif
    894944  if (quadraticObj) {
    895     // Use smaller
    896     actualDualStep_=min(actualDualStep_,actualPrimalStep_);
    897     actualPrimalStep_=actualDualStep_;
    898   }
     945    // Use smaller unless very small
     946    double smallerStep=min(actualDualStep_,actualPrimalStep_);
     947    if (smallerStep>0.0001) {
     948      actualDualStep_=smallerStep;
     949      actualPrimalStep_=smallerStep;
     950    }
     951  }
     952#define OFFQ
     953#ifndef OFFQ
     954  if (quadraticObj) {
     955    // Don't bother if phase 0 or 3 or large gap
     956    //if ((phase==1||phase==2||phase==0)&&maximumDualError_>0.1*complementarityGap_
     957    //&&smallerStep>0.001) {
     958    if ((phase==1||phase==2||phase==0||phase==3)) {
     959      // minimize complementarity + norm*dual inf ? primal inf
     960      // at first - just check better - if not
     961      // Complementarity gap will be a*change*change + b*change +c
     962      double a=0.0;
     963      double b=0.0;
     964      double c=0.0;
     965      /* SQUARE of dual infeasibility will be:
     966         square of dj - ......
     967      */
     968      double aq=0.0;
     969      double bq=0.0;
     970      double cq=0.0;
     971      double gamma2 = gamma_*gamma_; // gamma*gamma will be added to diagonal
     972      double * linearDjChange = new double[numberTotal];
     973      CoinZeroN(linearDjChange,numberColumns_);
     974      multiplyAdd(deltaY_,numberRows_,1.0,linearDjChange+numberColumns_,0.0);
     975      matrix_->transposeTimes(-1.0,deltaY_,linearDjChange);
     976      CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
     977      const int * columnQuadratic = quadratic->getIndices();
     978      const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
     979      const int * columnQuadraticLength = quadratic->getVectorLengths();
     980      double * quadraticElement = quadratic->getMutableElements();
     981      for (iColumn=0;iColumn<numberTotal;iColumn++) {
     982        double oldPrimal = solution_[iColumn];
     983        if (!flagged(iColumn)) {
     984          if (lowerBound(iColumn)) {
     985            double change =oldPrimal+deltaX_[iColumn]-lowerSlack_[iColumn]-lower_[iColumn];
     986            c += lowerSlack_[iColumn]*zVec_[iColumn];
     987            b += lowerSlack_[iColumn]*deltaZ_[iColumn]+zVec_[iColumn]*change;
     988            a += deltaZ_[iColumn]*change;
     989          }
     990          if (upperBound(iColumn)) {
     991            double change =upper_[iColumn]-oldPrimal-deltaX_[iColumn]-upperSlack_[iColumn];
     992            c += upperSlack_[iColumn]*wVec_[iColumn];
     993            b += upperSlack_[iColumn]*deltaW_[iColumn]+wVec_[iColumn]*change;
     994            a += deltaW_[iColumn]*change;
     995          }
     996          // new djs are dj_ + change*value
     997          double djChange = linearDjChange[iColumn];
     998          if (iColumn<numberColumns_) {
     999            for (CoinBigIndex j=columnQuadraticStart[iColumn];
     1000                 j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     1001              int jColumn = columnQuadratic[j];
     1002              double changeJ = deltaX_[jColumn];
     1003              double elementValue = quadraticElement[j];
     1004              djChange += changeJ*elementValue;
     1005            }
     1006          }
     1007          double gammaTerm = gamma2;
     1008          if (primalR_) {
     1009            gammaTerm += primalR_[iColumn];
     1010          }
     1011          djChange += gammaTerm;
     1012          // and dual infeasibility
     1013          double oldInf = dj_[iColumn]-zVec_[iColumn]+wVec_[iColumn]+
     1014            gammaTerm*solution_[iColumn];
     1015          double changeInf = djChange-deltaZ_[iColumn]+deltaW_[iColumn];
     1016          cq += oldInf*oldInf;
     1017          bq += 2.0*oldInf*changeInf;
     1018          aq += changeInf*changeInf;
     1019        } else {
     1020          // fixed
     1021          if (lowerBound(iColumn)) {
     1022            c += lowerSlack_[iColumn]*zVec_[iColumn];
     1023          }
     1024          if (upperBound(iColumn)) {
     1025            c += upperSlack_[iColumn]*wVec_[iColumn];
     1026          }
     1027          // new djs are dj_ + change*value
     1028          double djChange = linearDjChange[iColumn];
     1029          if (iColumn<numberColumns_) {
     1030            for (CoinBigIndex j=columnQuadraticStart[iColumn];
     1031                 j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     1032              int jColumn = columnQuadratic[j];
     1033              double changeJ = deltaX_[jColumn];
     1034              double elementValue = quadraticElement[j];
     1035              djChange += changeJ*elementValue;
     1036            }
     1037          }
     1038          double gammaTerm = gamma2;
     1039          if (primalR_) {
     1040            gammaTerm += primalR_[iColumn];
     1041          }
     1042          djChange += gammaTerm;
     1043          // and dual infeasibility
     1044          double oldInf = dj_[iColumn]-zVec_[iColumn]+wVec_[iColumn]+
     1045            gammaTerm*solution_[iColumn];
     1046          double changeInf = djChange-deltaZ_[iColumn]+deltaW_[iColumn];
     1047          cq += oldInf*oldInf;
     1048          bq += 2.0*oldInf*changeInf;
     1049          aq += changeInf*changeInf;
     1050        }
     1051      }
     1052      delete [] linearDjChange;
     1053      // ? We want to minimize complementarityGap + solutionNorm_*square of inf ??
     1054      // maybe use inf and do line search
     1055      // To check see if matches at current step
     1056      double step=actualPrimalStep_;
     1057      //Current gap + solutionNorm_ * sqrt (sum square inf)
     1058      double multiplier = solutionNorm_;
     1059      multiplier *= 0.01;
     1060      multiplier=1.0;
     1061      double currentInf =  multiplier*sqrt(cq);
     1062      double nextInf =  multiplier*sqrt(max(cq+step*bq+step*step*aq,0.0));
     1063      double allowedIncrease=1.4;
     1064#ifdef SOME_DEBUG
     1065      printf("lin %g %g %g -> %g\n",a,b,c,
     1066             c+b*step+a*step*step);
     1067      printf("quad %g %g %g -> %g\n",aq,bq,cq,
     1068             cq+bq*step+aq*step*step);
     1069      debugMove(7,step,step);
     1070      printf ("current dualInf %g, with step of %g is %g\n",
     1071              currentInf,step,nextInf);
     1072#endif
     1073      if (b>-1.0e-6) {
     1074        if (phase!=0)
     1075          directionNorm=-1.0;
     1076      }
     1077      if ((phase==1||phase==2||phase==0||phase==3)&&nextInf>0.1*complementarityGap_&&
     1078          nextInf>currentInf*allowedIncrease) {
     1079        //cq = max(cq,10.0);
     1080        // convert to (x+q)*(x+q) = w
     1081        double q = bq/(1.0*aq);
     1082        double w = max(q*q + (cq/aq)*(allowedIncrease-1.0),0.0);
     1083        w = sqrt(w);
     1084        double stepX = w-q;
     1085        step=stepX;
     1086        nextInf =
     1087          multiplier*sqrt(max(cq+step*bq+step*step*aq,0.0));
     1088#ifdef SOME_DEBUG
     1089        printf ("with step of %g dualInf is %g\n",
     1090                step,nextInf);
     1091#endif
     1092        actualDualStep_=min(step,actualDualStep_);
     1093        actualPrimalStep_=min(step,actualPrimalStep_);
     1094      }
     1095    }
     1096  } else {
     1097    // probably pointless as linear
     1098    // minimize complementarity
     1099    // Complementarity gap will be a*change*change + b*change +c
     1100    double a=0.0;
     1101    double b=0.0;
     1102    double c=0.0;
     1103    for (iColumn=0;iColumn<numberTotal;iColumn++) {
     1104      double oldPrimal = solution_[iColumn];
     1105      if (!flagged(iColumn)) {
     1106        if (lowerBound(iColumn)) {
     1107          double change =oldPrimal+deltaX_[iColumn]-lowerSlack_[iColumn]-lower_[iColumn];
     1108          c += lowerSlack_[iColumn]*zVec_[iColumn];
     1109          b += lowerSlack_[iColumn]*deltaZ_[iColumn]+zVec_[iColumn]*change;
     1110          a += deltaZ_[iColumn]*change;
     1111        }
     1112        if (upperBound(iColumn)) {
     1113          double change =upper_[iColumn]-oldPrimal-deltaX_[iColumn]-upperSlack_[iColumn];
     1114          c += upperSlack_[iColumn]*wVec_[iColumn];
     1115          b += upperSlack_[iColumn]*deltaW_[iColumn]+wVec_[iColumn]*change;
     1116          a += deltaW_[iColumn]*change;
     1117        }
     1118      } else {
     1119        // fixed
     1120        if (lowerBound(iColumn)) {
     1121          c += lowerSlack_[iColumn]*zVec_[iColumn];
     1122        }
     1123        if (upperBound(iColumn)) {
     1124          c += upperSlack_[iColumn]*wVec_[iColumn];
     1125        }
     1126      }
     1127    }
     1128    // ? We want to minimize complementarityGap;
     1129    // maybe use inf and do line search
     1130    // To check see if matches at current step
     1131    double step=min(actualPrimalStep_,actualDualStep_);
     1132    double next = c+b*step+a*step*step;
     1133#ifdef SOME_DEBUG
     1134    printf("lin %g %g %g -> %g\n",a,b,c,
     1135           c+b*step+a*step*step);
     1136    debugMove(7,step,step);
     1137#endif
     1138    if (b>-1.0e-6) {
     1139      if (phase==0) {
     1140#ifdef SOME_DEBUG
     1141        printf("*** odd phase 0 direction\n");
     1142#endif
     1143      } else {
     1144        directionNorm=-1.0;
     1145      }
     1146    }
     1147    // and with ratio
     1148    a=0.0;
     1149    b=0.0;
     1150    double ratio = actualDualStep_/actualPrimalStep_;
     1151    for (iColumn=0;iColumn<numberTotal;iColumn++) {
     1152      double oldPrimal = solution_[iColumn];
     1153      if (!flagged(iColumn)) {
     1154        if (lowerBound(iColumn)) {
     1155          double change =oldPrimal+deltaX_[iColumn]-lowerSlack_[iColumn]-lower_[iColumn];
     1156          b += lowerSlack_[iColumn]*deltaZ_[iColumn]*ratio+zVec_[iColumn]*change;
     1157          a += deltaZ_[iColumn]*change*ratio;
     1158        }
     1159        if (upperBound(iColumn)) {
     1160          double change =upper_[iColumn]-oldPrimal-deltaX_[iColumn]-upperSlack_[iColumn];
     1161          b += upperSlack_[iColumn]*deltaW_[iColumn]*ratio+wVec_[iColumn]*change;
     1162          a += deltaW_[iColumn]*change*ratio;
     1163        }
     1164      }
     1165    }
     1166    // ? We want to minimize complementarityGap;
     1167    // maybe use inf and do line search
     1168    // To check see if matches at current step
     1169    step=actualPrimalStep_;
     1170    double next2 = c+b*step+a*step*step;
     1171    if (next2>next) {
     1172      actualPrimalStep_=min(actualPrimalStep_,actualDualStep_);
     1173      actualDualStep_=actualPrimalStep_;
     1174    }
     1175#ifdef SOME_DEBUG
     1176    printf("linb %g %g %g -> %g\n",a,b,c,
     1177           c+b*step+a*step*step);
     1178    debugMove(7,actualPrimalStep_,actualDualStep_);
     1179#endif
     1180    if (b>-1.0e-6) {
     1181      if (phase==0) {
     1182#ifdef SOME_DEBUG
     1183        printf("*** odd phase 0 direction\n");
     1184#endif
     1185      } else {
     1186        directionNorm=-1.0;
     1187      }
     1188    }
     1189  }
     1190#else
     1191  //actualPrimalStep_ =0.5*actualDualStep_;
     1192#endif
    8991193#ifdef FULL_DEBUG
    9001194  if (phase==3){
     
    9051199        if (lowerBound(iColumn)) {
    9061200          double change = -rhsL_[iColumn] + deltaX_[iColumn];
    907           double dualValue=zVec[iColumn]+actualDualStep_*deltaZ_[iColumn];
    908           double primalValue=lowerSlack[iColumn]+actualPrimalStep_*change;
     1201          double dualValue=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn];
     1202          double primalValue=lowerSlack_[iColumn]+actualPrimalStep_*change;
    9091203          double gapProduct=dualValue*primalValue;
    9101204          if (delta2Z_[iColumn]<minBeta||delta2Z_[iColumn]>maxBeta)
     
    9141208        if (upperBound(iColumn)) {
    9151209          double change = rhsU_[iColumn]-deltaX_[iColumn];
    916           double dualValue=tVec[iColumn]+actualDualStep_*deltaT_[iColumn];
    917           double primalValue=upperSlack[iColumn]+actualPrimalStep_*change;
     1210          double dualValue=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn];
     1211          double primalValue=upperSlack_[iColumn]+actualPrimalStep_*change;
    9181212          double gapProduct=dualValue*primalValue;
    919           if (delta2T_[iColumn]<minBeta||delta2T_[iColumn]>maxBeta)
     1213          if (delta2W_[iColumn]<minBeta||delta2W_[iColumn]>maxBeta)
    9201214            printf("3upper %d primal %g, dual %g, gap %g, old gap %g\n",
    921                  iColumn,primalValue,dualValue,gapProduct,delta2T_[iColumn]);
     1215                 iColumn,primalValue,dualValue,gapProduct,delta2W_[iColumn]);
    9221216        }
    9231217      }
     
    9391233        if (lowerBound(iColumn)) {
    9401234          double change = -rhsL_[iColumn] + deltaX_[iColumn];
    941           double dualValue=zVec[iColumn]+actualDualStep_*deltaZ_[iColumn];
    942           double primalValue=lowerSlack[iColumn]+actualPrimalStep_*change;
     1235          double dualValue=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn];
     1236          double primalValue=lowerSlack_[iColumn]+actualPrimalStep_*change;
    9431237          double gapProduct=dualValue*primalValue;
    9441238          largestL = max(largestL,gapProduct);
     
    9491243        if (upperBound(iColumn)) {
    9501244          double change = rhsU_[iColumn]-deltaX_[iColumn];
    951           double dualValue=tVec[iColumn]+actualDualStep_*deltaT_[iColumn];
    952           double primalValue=upperSlack[iColumn]+actualPrimalStep_*change;
     1245          double dualValue=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn];
     1246          double primalValue=upperSlack_[iColumn]+actualPrimalStep_*change;
    9531247          double gapProduct=dualValue*primalValue;
    9541248          largestU = max(largestU,gapProduct);
     
    9711265        if (lowerBound(iColumn)) {
    9721266          double change = -rhsL_[iColumn] + deltaX_[iColumn];
    973           double dualValue=zVec[iColumn]+actualDualStep_*deltaZ_[iColumn];
    974           double primalValue=lowerSlack[iColumn]+actualPrimalStep_*change;
     1267          double dualValue=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn];
     1268          double primalValue=lowerSlack_[iColumn]+actualPrimalStep_*change;
    9751269          double gapProduct=dualValue*primalValue;
    9761270          if (gapProduct<minBeta)
     
    9841278        if (upperBound(iColumn)) {
    9851279          double change = rhsU_[iColumn]-deltaX_[iColumn];
    986           double dualValue=tVec[iColumn]+actualDualStep_*deltaT_[iColumn];
    987           double primalValue=upperSlack[iColumn]+actualPrimalStep_*change;
     1280          double dualValue=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn];
     1281          double primalValue=upperSlack_[iColumn]+actualPrimalStep_*change;
    9881282          double gapProduct=dualValue*primalValue;
    9891283          if (gapProduct<minBeta)
     
    10931387  }
    10941388  double * newError = new double [numberRows_];
    1095   double * workArray = workArray_;
    10961389  int numberTotal = numberRows_+numberColumns_;
    10971390  //if flagged then entries zero so can do
     
    11021395    int iColumn;
    11031396    for (iColumn=0;iColumn<numberTotal;iColumn++)
    1104       deltaX_[iColumn] = workArray[iColumn] - solution_[iColumn];
     1397      deltaX_[iColumn] = workArray_[iColumn] - solution_[iColumn];
    11051398    multiplyAdd(deltaX_+numberColumns_,numberRows_,-1.0,deltaY_,0.0);
    11061399    matrix_->times(1.0,deltaX_,deltaY_);
     
    11131406    int iColumn;
    11141407    for (iColumn=0;iColumn<numberTotal;iColumn++)
    1115       deltaX_[iColumn] = workArray[iColumn] - solution_[iColumn];
     1408      deltaX_[iColumn] = workArray_[iColumn] - solution_[iColumn];
    11161409    multiplyAdd(deltaX_+numberColumns_,numberRows_,-1.0,deltaY_,0.0);
    11171410    matrix_->times(1.0,deltaX_,deltaY_);
     
    11811474      for (iColumn=0;iColumn<numberTotal;iColumn++)
    11821475        deltaX_[iColumn] = deltaX_[iColumn]*diagonal_[iColumn]
    1183           -workArray[iColumn];
     1476          -workArray_[iColumn];
    11841477    } else {
    11851478      // KKT
    11861479      solveSystem(deltaX_, deltaY_,
    1187                   workArray,newError,region1Save,regionSave,lastError>1.0e-5);
     1480                  workArray_,newError,region1Save,regionSave,lastError>1.0e-5);
    11881481    }
    11891482    multiplyAdd(deltaX_+numberColumns_,numberRows_,-1.0,newError,0.0);
     
    12711564          for (iColumn=0;iColumn<numberTotal;iColumn++)
    12721565            deltaX_[iColumn] = deltaX_[iColumn]*diagonal_[iColumn]
    1273               -workArray[iColumn];
     1566              -workArray_[iColumn];
    12741567        } else {
    12751568          // KKT
     
    12881581  delete [] newError;
    12891582  // now rest
    1290   double * zVec = zVec_;
    1291   double * tVec = tVec_;
    1292   double * lowerSlack = lowerSlack_;
    1293   double * upperSlack = upperSlack_;
    12941583  double extra=eExtra;
    1295   //multiplyAdd(deltaY_,numberRows_,1.0,deltaT_+numberColumns_,0.0);
    1296   //CoinZeroN(deltaT_,numberColumns_);
    1297   //matrix_->transposeTimes(-1.0,deltaY_,deltaT_);
    1298 
     1584  //multiplyAdd(deltaY_,numberRows_,1.0,deltaW_+numberColumns_,0.0);
     1585  //CoinZeroN(deltaW_,numberColumns_);
     1586  //matrix_->transposeTimes(-1.0,deltaY_,deltaW_);
     1587 
    12991588  for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
    13001589    deltaSU_[iColumn]=0.0;
    13011590    deltaSL_[iColumn]=0.0;
    13021591    deltaZ_[iColumn]=0.0;
    1303     double dd=deltaT_[iColumn];
    1304     deltaT_[iColumn]=0.0;
     1592    double dd=deltaW_[iColumn];
     1593    deltaW_[iColumn]=0.0;
    13051594    if (!flagged(iColumn)) {
    13061595      double deltaX = deltaX_[iColumn];
    13071596      if (lowerBound(iColumn)) {
    13081597        double zValue = rhsZ_[iColumn];
    1309         double gHat = zValue + zVec[iColumn]*rhsL_[iColumn];
    1310         double slack = lowerSlack[iColumn]+extra;
     1598        double gHat = zValue + zVec_[iColumn]*rhsL_[iColumn];
     1599        double slack = lowerSlack_[iColumn]+extra;
    13111600        deltaSL_[iColumn] = -rhsL_[iColumn]+deltaX;
    1312         deltaZ_[iColumn]=(gHat-zVec[iColumn]*deltaX)/slack;
     1601        deltaZ_[iColumn]=(gHat-zVec_[iColumn]*deltaX)/slack;
    13131602      }
    13141603      if (upperBound(iColumn)) {
    1315         double tValue = rhsT_[iColumn];
    1316         double hHat = tValue - tVec[iColumn]*rhsU_[iColumn];
    1317         double slack = upperSlack[iColumn]+extra;
     1604        double wValue = rhsW_[iColumn];
     1605        double hHat = wValue - wVec_[iColumn]*rhsU_[iColumn];
     1606        double slack = upperSlack_[iColumn]+extra;
    13181607        deltaSU_[iColumn] = rhsU_[iColumn]-deltaX;
    1319         deltaT_[iColumn]=(hHat+tVec[iColumn]*deltaX)/slack;
     1608        deltaW_[iColumn]=(hHat+wVec_[iColumn]*deltaX)/slack;
    13201609      }
    13211610      if (0) {
     
    13231612        double gamma2 = gamma_*gamma_;
    13241613        double dZ=0.0;
    1325         double dT=0.0;
     1614        double dW=0.0;
    13261615        double zValue = rhsZ_[iColumn];
    1327         double gHat = zValue + zVec[iColumn]*rhsL_[iColumn];
    1328         double slackL = lowerSlack[iColumn]+extra;
    1329         double tValue = rhsT_[iColumn];
    1330         double hHat = tValue - tVec[iColumn]*rhsU_[iColumn];
    1331         double slackU = upperSlack[iColumn]+extra;
     1616        double gHat = zValue + zVec_[iColumn]*rhsL_[iColumn];
     1617        double slackL = lowerSlack_[iColumn]+extra;
     1618        double wValue = rhsW_[iColumn];
     1619        double hHat = wValue - wVec_[iColumn]*rhsU_[iColumn];
     1620        double slackU = upperSlack_[iColumn]+extra;
    13321621        double q = rhsC_[iColumn]+gamma2 * deltaX +dd;
    13331622        if (primalR_)
    13341623          q += deltaX*primalR_[iColumn];
    1335         dT = (gHat+hHat -slackL*q + (tValue-zValue)*deltaX)/(slackL+slackU);
    1336         dZ = dT + q;
     1624        dW = (gHat+hHat -slackL*q + (wValue-zValue)*deltaX)/(slackL+slackU);
     1625        dZ = dW + q;
    13371626        //printf("B %d old %g %g new %g %g\n",iColumn,deltaZ_[iColumn],
    1338         //deltaT_[iColumn],dZ,dT);
     1627        //deltaW_[iColumn],dZ,dW);
    13391628        if (lowerBound(iColumn)) {
    13401629          if (upperBound(iColumn)) {
    13411630            //printf("B %d old %g %g new %g %g\n",iColumn,deltaZ_[iColumn],
    1342             //deltaT_[iColumn],dZ,dT);
    1343             deltaT_[iColumn]=dT;
     1631            //deltaW_[iColumn],dZ,dW);
     1632            deltaW_[iColumn]=dW;
    13441633            deltaZ_[iColumn]=dZ;
    13451634          } else {
     
    13501639        } else {
    13511640          assert (upperBound(iColumn));
    1352           //printf("U %d old %g new %g\n",iColumn,deltaT_[iColumn],
    1353           //dT);
     1641          //printf("U %d old %g new %g\n",iColumn,deltaW_[iColumn],
     1642          //dW);
    13541643        }
    13551644      }
    13561645    }
    1357   }
     1646  }
     1647#if 0
     1648  double * check = new double[numberTotal]; 
     1649  // Check out rhsC_
     1650  multiplyAdd(deltaY_,numberRows_,-1.0,check+numberColumns_,0.0);
     1651  CoinZeroN(check,numberColumns_);
     1652  matrix_->transposeTimes(1.0,deltaY_,check);
     1653  quadraticDjs(check,deltaX_,-1.0);
     1654  for (iColumn=0;iColumn<numberTotal;iColumn++) {
     1655    check[iColumn] += deltaZ_[iColumn]-deltaW_[iColumn];
     1656    if (fabs(check[iColumn]-rhsC_[iColumn])>1.0e-3)
     1657      printf("rhsC %d %g %g\n",iColumn,check[iColumn],rhsC_[iColumn]);
     1658  }
     1659  // Check out rhsZ_
     1660  for (iColumn=0;iColumn<numberTotal;iColumn++) {
     1661    check[iColumn] += lowerSlack_[iColumn]*deltaZ_[iColumn]+
     1662      zVec_[iColumn]*deltaSL_[iColumn];
     1663    if (fabs(check[iColumn]-rhsZ_[iColumn])>1.0e-3)
     1664      printf("rhsZ %d %g %g\n",iColumn,check[iColumn],rhsZ_[iColumn]);
     1665  }
     1666  // Check out rhsW_
     1667  for (iColumn=0;iColumn<numberTotal;iColumn++) {
     1668    check[iColumn] += upperSlack_[iColumn]*deltaW_[iColumn]+
     1669      wVec_[iColumn]*deltaSU_[iColumn];
     1670    if (fabs(check[iColumn]-rhsW_[iColumn])>1.0e-3)
     1671      printf("rhsW %d %g %g\n",iColumn,check[iColumn],rhsW_[iColumn]);
     1672  }
     1673  delete [] check;
     1674#endif
    13581675  return relativeError;
    13591676}
     
    14211738      int numberColumns = quadratic->getNumCols();
    14221739      double scale = 1.0/scaleFactor;
    1423       if (scalingFlag_&&rowScale_) {
     1740      if (scalingFlag_>0&&rowScale_) {
    14241741        for (int iColumn=0;iColumn<numberColumns;iColumn++) {
    14251742          double scaleI = columnScale_[iColumn]*scale;
     
    15971914  }
    15981915  double gamma2 = gamma_*gamma_; // gamma*gamma will be added to diagonal
     1916#if 0
     1917  fakeSolution[0 ] =   0.072310129 ;
     1918  fakeSolution[1 ] =   0.053083871;
     1919  fakeSolution[2 ] =      0.178127;
     1920  fakeSolution[3 ] =    0.13215151;
     1921  fakeSolution[4 ] =   0.072715642;
     1922  fakeSolution[5 ] =    0.15680727;
     1923  fakeSolution[6 ] =    0.16841689;
     1924  fakeSolution[7 ] =   0.093612798 ;
     1925  fakeSolution[8 ] =   0.072774891 ;
     1926  fakeSolution[9]=1.0;
     1927  initialValue=1.0e-5;
     1928  safeObjectiveValue=1.0e-5;
     1929#endif
    15991930  // First do primal side
    16001931  for ( iColumn=0;iColumn<numberTotal;iColumn++) {
     
    16601991      solution_[iColumn]=lower_[iColumn];
    16611992      zVec_[iColumn]=0.0;
    1662       tVec_[iColumn]=0.0;
     1993      wVec_[iColumn]=0.0;
    16631994      diagonal_[iColumn]=0.0;
    16641995    }
     
    17222053      double upperValue=upper_[iColumn];
    17232054      double reducedCost=dj_[iColumn];
    1724       reducedCost = cost_[iColumn];
    17252055      double low=0.0;
    17262056      double high=0.0;
     
    17432073          }
    17442074          double t = high+extra;
    1745           double ratioW;
     2075          double ratioT;
    17462076          if (t<zwLarge) {
    1747             ratioW=1.0;
     2077            ratioT=1.0;
    17482078          } else {
    1749             ratioW=sqrt(zwLarge/t);
     2079            ratioT=sqrt(zwLarge/t);
    17502080          }
    17512081          //modify s and t
     
    17592089          if (reducedCost>=0.0) {
    17602090            zVec_[iColumn]=reducedCost + safeObjectiveValue*ratioZ;
    1761             tVec_[iColumn]=safeObjectiveValue*ratioW;
     2091            zVec_[iColumn]=max(reducedCost, safeObjectiveValue*ratioZ);
     2092            wVec_[iColumn]=safeObjectiveValue*ratioT;
    17622093          } else {
    17632094            zVec_[iColumn]=safeObjectiveValue*ratioZ;
    1764             tVec_[iColumn]=-reducedCost + safeObjectiveValue*ratioW;
     2095            wVec_[iColumn]=-reducedCost + safeObjectiveValue*ratioT;
     2096            wVec_[iColumn]=max(-reducedCost , safeObjectiveValue*ratioT);
    17652097          }
    17662098          double gammaTerm = gamma2;
     
    17682100            gammaTerm += primalR_[iColumn];
    17692101          diagonal_[iColumn] = (t*s)/
    1770             (s*tVec_[iColumn]+t*zVec_[iColumn]+gammaTerm*t*s);
     2102            (s*wVec_[iColumn]+t*zVec_[iColumn]+gammaTerm*t*s);
    17712103        } else {
    17722104          //just lower bound
     
    17862118          if (reducedCost>=0.0) {
    17872119            zVec_[iColumn]=reducedCost + safeObjectiveValue*ratioZ;
    1788             tVec_[iColumn]=0.0;
     2120            zVec_[iColumn]=max(reducedCost , safeObjectiveValue*ratioZ);
     2121            wVec_[iColumn]=0.0;
    17892122          } else {
    17902123            zVec_[iColumn]=safeObjectiveValue*ratioZ;
    1791             tVec_[iColumn]=0.0;
     2124            wVec_[iColumn]=0.0;
    17922125          }
    17932126          double gammaTerm = gamma2;
     
    18022135          high=upperValue-primalValue;
    18032136          double t = high+extra;
    1804           double ratioW;
     2137          double ratioT;
    18052138          if (t<zwLarge) {
    1806             ratioW=1.0;
     2139            ratioT=1.0;
    18072140          } else {
    1808             ratioW=sqrt(zwLarge/t);
     2141            ratioT=sqrt(zwLarge/t);
    18092142          }
    18102143          //modify t
     
    18142147          if (reducedCost>=0.0) {
    18152148            zVec_[iColumn]=0.0;
    1816             tVec_[iColumn]=safeObjectiveValue*ratioW;
     2149            wVec_[iColumn]=safeObjectiveValue*ratioT;
    18172150          } else {
    18182151            zVec_[iColumn]=0.0;
    1819             tVec_[iColumn]=-reducedCost + safeObjectiveValue*ratioW;
     2152            wVec_[iColumn]=-reducedCost + safeObjectiveValue*ratioT;
     2153            wVec_[iColumn]=max(-reducedCost , safeObjectiveValue*ratioT);
    18202154          }
    18212155          double gammaTerm = gamma2;
    18222156          if (primalR_)
    18232157            gammaTerm += primalR_[iColumn];
    1824           diagonal_[iColumn] =  t/(tVec_[iColumn]+t*gammaTerm);
     2158          diagonal_[iColumn] =  t/(wVec_[iColumn]+t*gammaTerm);
    18252159        }
    18262160      }
     
    18352169             diagonal_[i],fabs(dj_[i]),
    18362170             lowerSlack_[i],zVec_[i],
    1837              upperSlack_[i],tVec_[i]);
     2171             upperSlack_[i],wVec_[i]);
    18382172  } else {
    18392173    for (int i=0;i<numberTotal;i++)
    18402174      printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n",i,fabs(solution_[i]),
    18412175             diagonal_[i],fabs(dj_[i]),
    1842              upperSlack_[i],tVec_[i],
     2176             upperSlack_[i],wVec_[i],
    18432177             lowerSlack_[i],zVec_[i] );
    18442178  }
     
    18722206  double primalTolerance =  dblParam_[ClpPrimalTolerance];
    18732207  dualTolerance=dualTolerance/scaleFactor_;
    1874   double * zVec = zVec_;
    1875   double * tVec = tVec_;
    1876   double * primal = solution_;
    1877   double * lower = lower_;
    1878   double * upper = upper_;
    1879   double * lowerSlack = lowerSlack_;
    1880   double * upperSlack = upperSlack_;
    1881   double * deltaZ = deltaZ_;
    1882   double * deltaT = deltaT_;
    1883   double * deltaX = deltaX_;
    18842208  for (int iColumn=0;iColumn<numberTotal;iColumn++) {
    18852209    if (!fixedOrFree(iColumn)) {
     
    18872211      //can collapse as if no lower bound both zVec and deltaZ 0.0
    18882212      double newZ=0.0;
    1889       double newT=0.0;
     2213      double newW=0.0;
    18902214      if (lowerBound(iColumn)) {
    18912215        numberComplementarityItems++;
     
    18932217        double primalValue;
    18942218        if (!phase) {
    1895           dualValue=zVec[iColumn];
    1896           primalValue=lowerSlack[iColumn];
     2219          dualValue=zVec_[iColumn];
     2220          primalValue=lowerSlack_[iColumn];
    18972221        } else {
    18982222          double change;
    1899           change =primal[iColumn]+deltaX[iColumn]-lowerSlack[iColumn]-lower[iColumn];
    1900           dualValue=zVec[iColumn]+actualDualStep_*deltaZ[iColumn];
     2223          change =solution_[iColumn]+deltaX_[iColumn]-lowerSlack_[iColumn]-lower_[iColumn];
     2224          dualValue=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn];
    19012225          newZ=dualValue;
    1902           primalValue=lowerSlack[iColumn]+actualPrimalStep_*change;
    1903         } 
     2226          primalValue=lowerSlack_[iColumn]+actualPrimalStep_*change;
     2227        }
    19042228        //reduce primalValue
    19052229        if (primalValue>largeGap) {
     
    19152239        }
    19162240        gap+=gapProduct;
     2241        //printf("l %d prim %g dual %g totalGap %g\n",
     2242        //   iColumn,primalValue,dualValue,gap);
    19172243        if (gapProduct>largestGap) {
    19182244          largestGap=gapProduct;
     
    19282254        double primalValue;
    19292255        if (!phase) {
    1930           dualValue=tVec[iColumn];
    1931           primalValue=upperSlack[iColumn];
     2256          dualValue=wVec_[iColumn];
     2257          primalValue=upperSlack_[iColumn];
    19322258        } else {
    19332259          double change;
    1934           change =upper[iColumn]-primal[iColumn]-deltaX[iColumn]-upperSlack[iColumn];
    1935           dualValue=tVec[iColumn]+actualDualStep_*deltaT[iColumn];
    1936           newT=dualValue;
    1937           primalValue=upperSlack[iColumn]+actualPrimalStep_*change;
     2260          change =upper_[iColumn]-solution_[iColumn]-deltaX_[iColumn]-upperSlack_[iColumn];
     2261          dualValue=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn];
     2262          newW=dualValue;
     2263          primalValue=upperSlack_[iColumn]+actualPrimalStep_*change;
    19382264        }
    19392265        //reduce primalValue
     
    19502276        }
    19512277        gap+=gapProduct;
     2278        //printf("u %d prim %g dual %g totalGap %g\n",
     2279        //   iColumn,primalValue,dualValue,gap);
    19522280        if (gapProduct>largestGap) {
    19532281          largestGap=gapProduct;
     
    19592287    }
    19602288  }
     2289  //if (numberIterations_>4)
     2290  //exit(9);
    19612291  if (!phase&&numberNegativeGaps) {
    19622292      handler_->message(CLP_BARRIER_NEGATIVE_GAPS,messages_)
     
    19692299    numberComplementarityPairs=1;
    19702300  }
    1971   //printf("gap %g - smallest %g, largest %g, pairs %d\n",
    1972   // gap,smallestGap,largestGap,numberComplementarityPairs);
     2301#ifdef SOME_DEBUG
     2302  printf("with d,p steps %g,%g gap %g - smallest %g, largest %g, pairs %d\n",
     2303         actualDualStep_,actualPrimalStep_,
     2304         gap,smallestGap,largestGap,numberComplementarityPairs);
     2305#endif
    19732306  return gap;
    19742307}
     
    19782311{
    19792312  double extra =eExtra;
    1980   double * zVec = zVec_;
    1981   double * tVec = tVec_;
    1982   double * primal = solution_;
    1983   double * dj = dj_;
    1984   double * lower = lower_;
    1985   double * upper = upper_;
    1986   double * lowerSlack = lowerSlack_;
    1987   double * upperSlack = upperSlack_;
    1988   double * workArray = workArray_;
    19892313  int numberTotal = numberRows_ + numberColumns_;
    19902314  int iColumn;
     
    20102334      rhsL_[iColumn]=0.0;
    20112335      rhsZ_[iColumn]=0.0;
    2012       rhsT_[iColumn]=0.0;
     2336      rhsW_[iColumn]=0.0;
    20132337      if (!flagged(iColumn)) {
    2014         rhsC_[iColumn] = dj[iColumn]-zVec[iColumn]+tVec[iColumn];
     2338        rhsC_[iColumn] = dj_[iColumn]-zVec_[iColumn]+wVec_[iColumn];
    20152339        rhsC_[iColumn] += gamma2*solution_[iColumn];
    20162340        if (primalR_)
    20172341          rhsC_[iColumn] += primalR_[iColumn]*solution_[iColumn];
    20182342        if (lowerBound(iColumn)) {
    2019           rhsZ_[iColumn] = -zVec[iColumn]*(lowerSlack[iColumn]+extra);
    2020           rhsL_[iColumn] = max(0.0,(lower[iColumn]+lowerSlack[iColumn])-solution_[iColumn]);
     2343          rhsZ_[iColumn] = -zVec_[iColumn]*(lowerSlack_[iColumn]+extra);
     2344          rhsL_[iColumn] = max(0.0,(lower_[iColumn]+lowerSlack_[iColumn])-solution_[iColumn]);
    20212345        }
    20222346        if (upperBound(iColumn)) {
    2023           rhsT_[iColumn] = -tVec[iColumn]*(upperSlack[iColumn]+extra);
    2024           rhsU_[iColumn] = min(0.0,(upper[iColumn]-upperSlack[iColumn])-primal[iColumn]);
     2347          rhsW_[iColumn] = -wVec_[iColumn]*(upperSlack_[iColumn]+extra);
     2348          rhsU_[iColumn] = min(0.0,(upper_[iColumn]-upperSlack_[iColumn])-solution_[iColumn]);
    20252349        }
    20262350      }
     
    20312355             diagonal_[i],dj_[i],
    20322356             lowerSlack_[i],zVec_[i],
    2033              upperSlack_[i],tVec_[i]);
     2357             upperSlack_[i],wVec_[i]);
    20342358      printf("%d %.18g %.18g %.18g %.18g %.18g\n",i,rhsC_[i],
    20352359             rhsZ_[i],rhsL_[i],
    2036              rhsT_[i],rhsU_[i]);
     2360             rhsW_[i],rhsU_[i]);
    20372361    }
    20382362#if 0
     
    20402364      if (!fabs(rhsZ_[i]))
    20412365        rhsZ_[i]=0.0;
    2042       if (!fabs(rhsT_[i]))
    2043         rhsT_[i]=0.0;
     2366      if (!fabs(rhsW_[i]))
     2367        rhsW_[i]=0.0;
    20442368      if (!fabs(rhsU_[i]))
    20452369        rhsU_[i]=0.0;
     
    20522376               diagonal_[i],dj_[i],
    20532377               lowerSlack_[i],zVec_[i],
    2054                upperSlack_[i],tVec_[i]);
     2378               upperSlack_[i],wVec_[i]);
    20552379      for (int i=0;i<3;i++)
    20562380        printf("%d %.18g %.18g %.18g %.18g %.18g\n",i,rhsC_[i],
    20572381               rhsZ_[i],rhsL_[i],
    2058                rhsT_[i],rhsU_[i]);
     2382               rhsW_[i],rhsU_[i]);
    20592383    } else {
    20602384      for (int i=0;i<3;i++)
     
    20622386               diagonal_[i],dj_[i],
    20632387               lowerSlack_[i],zVec_[i],
    2064                upperSlack_[i],tVec_[i]);
     2388               upperSlack_[i],wVec_[i]);
    20652389      for (int i=0;i<3;i++)
    20662390        printf("%d %.18g %.18g %.18g %.18g %.18g\n",i,rhsC_[i],
    20672391               rhsZ_[i],rhsL_[i],
    2068                rhsT_[i],rhsU_[i]);
     2392               rhsW_[i],rhsU_[i]);
    20692393    }
    20702394#endif
     
    20742398    for (iColumn=0;iColumn<numberTotal;iColumn++) {
    20752399      rhsZ_[iColumn]=0.0;
    2076       rhsT_[iColumn]=0.0;
     2400      rhsW_[iColumn]=0.0;
    20772401      if (!flagged(iColumn)) {
    20782402        if (lowerBound(iColumn)) {
    2079           rhsZ_[iColumn] = mu_ -zVec[iColumn]*(lowerSlack[iColumn]+extra)
     2403          rhsZ_[iColumn] = mu_ -zVec_[iColumn]*(lowerSlack_[iColumn]+extra)
    20802404            - deltaZ_[iColumn]*deltaX_[iColumn];
    20812405          // To bring in line with OSL
     
    20832407        }
    20842408        if (upperBound(iColumn)) {
    2085           rhsT_[iColumn] = mu_ -tVec[iColumn]*(upperSlack[iColumn]+extra)
    2086             +deltaT_[iColumn]*deltaX_[iColumn];
     2409          rhsW_[iColumn] = mu_ -wVec_[iColumn]*(upperSlack_[iColumn]+extra)
     2410            +deltaW_[iColumn]*deltaX_[iColumn];
    20872411          // To bring in line with OSL
    2088           rhsT_[iColumn] -= deltaT_[iColumn]*rhsU_[iColumn];
     2412          rhsW_[iColumn] -= deltaW_[iColumn]*rhsU_[iColumn];
    20892413        }
    20902414      }
     
    20942418      if (!fabs(rhsZ_[i]))
    20952419        rhsZ_[i]=0.0;
    2096       if (!fabs(rhsT_[i]))
    2097         rhsT_[i]=0.0;
     2420      if (!fabs(rhsW_[i]))
     2421        rhsW_[i]=0.0;
    20982422      if (!fabs(rhsU_[i]))
    20992423        rhsU_[i]=0.0;
     
    21062430               diagonal_[i],fabs(dj_[i]),
    21072431               lowerSlack_[i],zVec_[i],
    2108                upperSlack_[i],tVec_[i]);
     2432               upperSlack_[i],wVec_[i]);
    21092433      for (int i=0;i<numberTotal;i++)
    21102434        printf("%d %.18g %.18g %.18g %.18g %.18g\n",i,fabs(rhsC_[i]),
    21112435               rhsZ_[i],rhsL_[i],
    2112                rhsT_[i],rhsU_[i]);
     2436               rhsW_[i],rhsU_[i]);
    21132437    } else {
    21142438      for (int i=0;i<numberTotal;i++)
    21152439        printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n",i,fabs(solution_[i]),
    21162440               diagonal_[i],fabs(dj_[i]),
    2117                upperSlack_[i],tVec_[i],
     2441               upperSlack_[i],wVec_[i],
    21182442               lowerSlack_[i],zVec_[i] );
    21192443      for (int i=0;i<numberTotal;i++)
    21202444        printf("%d %.18g %.18g %.18g %.18g %.18g\n",i,fabs(rhsC_[i]),
    2121                rhsT_[i],rhsU_[i],
     2445               rhsW_[i],rhsU_[i],
    21222446               rhsZ_[i],rhsL_[i]);
    21232447    }
     
    21292453    for (iColumn=0;iColumn<numberTotal;iColumn++) {
    21302454      rhsZ_[iColumn]=0.0;
    2131       rhsT_[iColumn]=0.0;
     2455      rhsW_[iColumn]=0.0;
    21322456      if (!flagged(iColumn)) {
    21332457        if (lowerBound(iColumn)) {
    2134           rhsZ_[iColumn] = mu_ - zVec[iColumn]*(lowerSlack[iColumn]+extra);
     2458          rhsZ_[iColumn] = mu_ - zVec_[iColumn]*(lowerSlack_[iColumn]+extra);
    21352459        }
    21362460        if (upperBound(iColumn)) {
    2137           rhsT_[iColumn] = mu_ - tVec[iColumn]*(upperSlack[iColumn]+extra);
     2461          rhsW_[iColumn] = mu_ - wVec_[iColumn]*(upperSlack_[iColumn]+extra);
    21382462        }
    21392463      }
     
    21552479          if (lowerBound(iColumn)) {
    21562480            double change = -rhsL_[iColumn] + deltaX_[iColumn];
    2157             double dualValue=zVec[iColumn]+dualStep*deltaZ_[iColumn];
    2158             double primalValue=lowerSlack[iColumn]+primalStep*change;
     2481            double dualValue=zVec_[iColumn]+dualStep*deltaZ_[iColumn];
     2482            double primalValue=lowerSlack_[iColumn]+primalStep*change;
    21592483            double gapProduct=dualValue*primalValue;
    21602484            if (gapProduct>0.0&&dualValue<0.0)
     
    21802504          if (upperBound(iColumn)) {
    21812505            double change = rhsU_[iColumn]-deltaX_[iColumn];
    2182             double dualValue=tVec[iColumn]+dualStep*deltaT_[iColumn];
    2183             double primalValue=upperSlack[iColumn]+primalStep*change;
     2506            double dualValue=wVec_[iColumn]+dualStep*deltaW_[iColumn];
     2507            double primalValue=upperSlack_[iColumn]+primalStep*change;
    21842508            double gapProduct=dualValue*primalValue;
    21852509            if (gapProduct>0.0&&dualValue<0.0)
    21862510              gapProduct = - gapProduct;
    21872511#ifdef FULL_DEBUG
    2188             delta2T_[iColumn]=gapProduct;
    2189             if (delta2T_[iColumn]<minBeta||delta2T_[iColumn]>maxBeta)
     2512            delta2W_[iColumn]=gapProduct;
     2513            if (delta2W_[iColumn]<minBeta||delta2W_[iColumn]>maxBeta)
    21902514              printf("upper %d primal %g, dual %g, gap %g\n",
    21912515                     iColumn,primalValue,dualValue,gapProduct);
     
    21992523              assert (value<0.0);
    22002524            }
    2201             rhsT_[iColumn] += value;
     2525            rhsW_[iColumn] += value;
    22022526          }
    22032527        }
     
    22102534      double value = rhsC_[iColumn];
    22112535      double zValue = rhsZ_[iColumn];
    2212       double tValue = rhsT_[iColumn];
     2536      double wValue = rhsW_[iColumn];
    22132537#if 0
    22142538#if 1
     
    22172541        value = dj[iColumn];
    22182542        zValue=0.0;
    2219         tValue=0.0;
     2543        wValue=0.0;
    22202544      } else if (phase==2) {
    22212545        // more accurate
    22222546        value = dj[iColumn];
    22232547        zValue=mu_;
    2224         tValue=mu_;
     2548        wValue=mu_;
    22252549      }
    22262550#endif
     
    22282552      assert (rhsU_[iColumn]<=0.0);
    22292553      if (lowerBound(iColumn)) {
    2230         value += (-zVec[iColumn]*rhsL_[iColumn]-zValue)/
    2231           (lowerSlack[iColumn]+extra);
     2554        value += (-zVec_[iColumn]*rhsL_[iColumn]-zValue)/
     2555          (lowerSlack_[iColumn]+extra);
    22322556      }
    22332557      if (upperBound(iColumn)) {
    2234         value += (tValue-tVec[iColumn]*rhsU_[iColumn])/
    2235           (upperSlack[iColumn]+extra);
     2558        value += (wValue-wVec_[iColumn]*rhsU_[iColumn])/
     2559          (upperSlack_[iColumn]+extra);
    22362560      }
    22372561#else
    22382562      if (lowerBound(iColumn)) {
    2239         double gHat = zValue + zVec[iColumn]*rhsL_[iColumn];
    2240         value -= gHat/(lowerSlack[iColumn]+extra);
     2563        double gHat = zValue + zVec_[iColumn]*rhsL_[iColumn];
     2564        value -= gHat/(lowerSlack_[iColumn]+extra);
    22412565      }
    22422566      if (upperBound(iColumn)) {
    2243         double hHat = tValue - tVec[iColumn]*rhsU_[iColumn];
    2244         value += hHat/(upperSlack[iColumn]+extra);
    2245       }
    2246 #endif
    2247       workArray[iColumn]=diagonal_[iColumn]*value;
     2567        double hHat = wValue - wVec_[iColumn]*rhsU_[iColumn];
     2568        value += hHat/(upperSlack_[iColumn]+extra);
     2569      }
     2570#endif
     2571      workArray_[iColumn]=diagonal_[iColumn]*value;
    22482572    }
    22492573#if 0
    22502574    if (solution_[0]>0.0) {
    22512575      for (int i=0;i<numberTotal;i++)
    2252         printf("%d %.18g\n",i,workArray[i]);
     2576        printf("%d %.18g\n",i,workArray_[i]);
    22532577    } else {
    22542578      for (int i=0;i<numberTotal;i++)
    2255         printf("%d %.18g\n",i,workArray[i]);
     2579        printf("%d %.18g\n",i,workArray_[i]);
    22562580    }
    22572581    exit(66);
     
    22622586      double value = rhsC_[iColumn];
    22632587      double zValue = rhsZ_[iColumn];
    2264       double tValue = rhsT_[iColumn];
     2588      double wValue = rhsW_[iColumn];
    22652589      if (lowerBound(iColumn)) {
    2266         double gHat = zValue + zVec[iColumn]*rhsL_[iColumn];
    2267         value -= gHat/(lowerSlack[iColumn]+extra);
     2590        double gHat = zValue + zVec_[iColumn]*rhsL_[iColumn];
     2591        value -= gHat/(lowerSlack_[iColumn]+extra);
    22682592      }
    22692593      if (upperBound(iColumn)) {
    2270         double hHat = tValue - tVec[iColumn]*rhsU_[iColumn];
    2271         value += hHat/(upperSlack[iColumn]+extra);
    2272       }
    2273       workArray[iColumn]=value;
     2594        double hHat = wValue - wVec_[iColumn]*rhsU_[iColumn];
     2595        value += hHat/(upperSlack_[iColumn]+extra);
     2596      }
     2597      workArray_[iColumn]=value;
    22742598    }
    22752599  }
     
    23312655  if (goodMove)
    23322656    goodMove=checkGoodMove2(step,bestNextGap,allowIncreasingGap);
     2657  // Say good if small
     2658  //if (quadraticObj) {
     2659  if (max(actualDualStep_,actualPrimalStep_)<1.0e-6)
     2660    goodMove=true;
    23332661  if (!goodMove) {
    23342662    //try smaller of two
     
    23422670    }
    23432671    actualPrimalStep_=step;
     2672    //if (quadraticObj)
     2673    //actualPrimalStep_ *=0.5;
    23442674    actualDualStep_=step;
    23452675    goodMove=checkGoodMove2(step,bestNextGap,allowIncreasingGap);
     2676    int pass=0;
    23462677    while (!goodMove) {
     2678      pass++;
    23472679      double gap = bestNextGap;
    23482680      goodMove=checkGoodMove2(step,gap,allowIncreasingGap);
    2349       if (goodMove) {
     2681      if (goodMove||pass>3) {
    23502682        returnGap=gap;
    23512683        break;
    23522684      }
    2353       if (step<1.0e-10) {
     2685      if (step<1.0e-4) {
    23542686        break;
    23552687      }
    23562688      step*=0.5;
    23572689      actualPrimalStep_=step;
     2690      //if (quadraticObj)
     2691      //actualPrimalStep_ *=0.5;
    23582692      actualDualStep_=step;
    23592693    } /* endwhile */
     
    23852719    CoinMemcpyN(deltaY_,numberRows_,workArray+numberColumns_);
    23862720    matrix_->transposeTimes(-1.0,deltaY_,workArray);
    2387     double * deltaZ = deltaZ_;
    2388     double * deltaT = deltaT_;
    2389     double * lower = lower_;
    2390     double * upper = upper_;
    2391     //direction vector in deltaX
    2392     double * deltaX = deltaX_;
    2393     double * cost = cost_;
    23942721    //double sumPerturbCost=0.0;
    23952722    for (int iColumn=0;iColumn<numberTotal;iColumn++) {
    23962723      if (!flagged(iColumn)) {
    23972724        if (lowerBound(iColumn)) {
    2398           //sumPerturbCost+=deltaX[iColumn];
    2399           deltaObjectiveDual+=deltaZ[iColumn]*lower[iColumn];
     2725          //sumPerturbCost+=deltaX_[iColumn];
     2726          deltaObjectiveDual+=deltaZ_[iColumn]*lower_[iColumn];
    24002727        }
    24012728        if (upperBound(iColumn)) {
    2402           //sumPerturbCost-=deltaX[iColumn];
    2403           deltaObjectiveDual-=deltaT[iColumn]*upper[iColumn];
     2729          //sumPerturbCost-=deltaX_[iColumn];
     2730          deltaObjectiveDual-=deltaW_[iColumn]*upper_[iColumn];
    24042731        }
    2405         double change = fabs(workArray[iColumn]-deltaZ[iColumn]+deltaT[iColumn]);
     2732        double change = fabs(workArray_[iColumn]-deltaZ_[iColumn]+deltaW_[iColumn]);
    24062733        error = max (change,error);
    24072734      }
    2408       deltaObjectivePrimal += cost[iColumn] * deltaX[iColumn];
     2735      deltaObjectivePrimal += cost_[iColumn] * deltaX_[iColumn];
    24092736    }
    24102737    //deltaObjectivePrimal+=sumPerturbCost*linearPerturbation_;
    24112738    double testValue;
    24122739    if (error>0.0) {
    2413       testValue=1.0e1*maximumDualError_/error;
     2740      testValue=1.0e1*max(maximumDualError_,1.0e-12)/error;
    24142741    } else {
    24152742      testValue=1.0e1;
    24162743    }
    2417     if (testValue<actualDualStep_) {
     2744    // If quadratic then primal step may compensate
     2745    if (testValue<actualDualStep_&&!quadraticObj) {
    24182746      handler_->message(CLP_BARRIER_REDUCING,messages_)
    24192747      <<"dual"<<actualDualStep_
     
    24272755    //check change in AX not too much
    24282756    //??? could be dropped row going infeasible
    2429     double ratio = 1.0e1*maximumRHSError_/maximumRHSChange_;
     2757    double ratio = 1.0e1*max(maximumRHSError_,1.0e-12)/maximumRHSChange_;
    24302758    if (ratio<actualPrimalStep_) {
    24312759      handler_->message(CLP_BARRIER_REDUCING,messages_)
     
    24462774}
    24472775//:  checks for one step size
    2448 bool ClpPredictorCorrector::checkGoodMove2(const double move,
     2776bool ClpPredictorCorrector::checkGoodMove2(double move,
    24492777                                           double & bestNextGap,
    24502778                                           bool allowIncreasingGap)
     
    24532781  const double gamma = 1.0e-8;
    24542782  const double gammap = 1.0e-8;
    2455   const double gammad = 1.0e-8;
     2783  double gammad = 1.0e-8;
    24562784  int nextNumber;
    24572785  int nextNumberItems;
     
    24612789  double lowerBoundGap = gamma*nextGap*complementarityMultiplier;
    24622790  bool goodMove=true;
    2463   double * deltaZ = deltaZ_;
    2464   double * deltaT = deltaT_;
    2465   double * deltaSL = deltaSL_;
    2466   double * deltaSU = deltaSU_;
    2467   double * zVec = zVec_;
    2468   double * tVec = tVec_;
    2469   double * lowerSlack = lowerSlack_;
    2470   double * upperSlack = upperSlack_;
    24712791  for (int iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
    24722792    if (!flagged(iColumn)) {
    24732793      if (lowerBound(iColumn)) {
    2474         double part1=lowerSlack[iColumn]+actualPrimalStep_*deltaSL[iColumn];
    2475         double part2=zVec[iColumn]+actualDualStep_*deltaZ[iColumn];
     2794        double part1=lowerSlack_[iColumn]+actualPrimalStep_*deltaSL_[iColumn];
     2795        double part2=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn];
    24762796        if (part1*part2<lowerBoundGap) {
    24772797          goodMove=false;
     
    24802800      }
    24812801      if (upperBound(iColumn)) {
    2482         double part1=upperSlack[iColumn]+actualPrimalStep_*deltaSU[iColumn];
    2483         double part2=tVec[iColumn]+actualDualStep_*deltaT[iColumn];
     2802        double part1=upperSlack_[iColumn]+actualPrimalStep_*deltaSU_[iColumn];
     2803        double part2=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn];
    24842804        if (part1*part2<lowerBoundGap) {
    24852805          goodMove=false;
     
    24902810  }
    24912811   double * nextDj=NULL;
     2812   double maximumDualError = maximumDualError_;
    24922813#ifndef NO_RTTI
    24932814  ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
     
    24982819#endif
    24992820  if (quadraticObj) {
     2821    // change gammad
     2822    gammad=1.0e-4;
    25002823    double gamma2 = gamma_*gamma_;
    25012824    nextDj = new double [numberColumns_];
     
    25192842      if (!fixedOrFree(iColumn)) {
    25202843        double newZ=0.0;
    2521         double newT=0.0;
     2844        double newW=0.0;
    25222845        if (lowerBound(iColumn)) {
    2523           newZ=zVec[iColumn]+actualDualStep_*deltaZ[iColumn];
     2846          newZ=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn];
    25242847        }
    25252848        if (upperBound(iColumn)) {
    2526           newT=tVec[iColumn]+actualDualStep_*deltaT[iColumn];
     2849          newW=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn];
    25272850        }
    25282851        if (columnQuadraticLength[iColumn]) {
     
    25302853          if (primalR_)
    25312854            gammaTerm += primalR_[iColumn];
    2532           double dualInfeasibility=
    2533             dj_[iColumn]-zVec[iColumn]+tVec[iColumn]
    2534             +gammaTerm*solution_[iColumn];
     2855          //double dualInfeasibility=
     2856          //dj_[iColumn]-zVec_[iColumn]+wVec_[iColumn]
     2857          //+gammaTerm*solution_[iColumn];
    25352858          double newInfeasibility=
    2536             nextDj[iColumn]-newZ+newT
     2859            nextDj[iColumn]-newZ+newW
    25372860            +gammaTerm*(solution_[iColumn]+actualPrimalStep_*deltaX_[iColumn]);
    2538           if (fabs(newInfeasibility)>max(2000.0*maximumDualError_,1.0e-2)) {
    2539             if (dualInfeasibility*newInfeasibility<0.0) {
    2540               //printf("%d current %g next %g\n",iColumn,dualInfeasibility,
    2541               //     newInfeasibility);
    2542               goodMove=false;
    2543             }
    2544           }
     2861          maximumDualError = max(maximumDualError,newInfeasibility);
     2862          //if (fabs(newInfeasibility)>max(2000.0*maximumDualError_,1.0e-2)) {
     2863          //if (dualInfeasibility*newInfeasibility<0.0) {
     2864          //  printf("%d current %g next %g\n",iColumn,dualInfeasibility,
     2865          //       newInfeasibility);
     2866          //  goodMove=false;
     2867          //}
     2868          //}
    25452869        }
    25462870      }
     
    25562880    errorCheck=maximumBoundInfeasibility_;
    25572881  }
     2882  // scale back move
     2883  move = min(move,0.95);
    25582884  //scale
    25592885  if ((1.0-move)*errorCheck>primalTolerance()) {
     
    25632889  }
    25642890  //      Satisfy g_d(alpha)?
    2565   errorCheck=maximumDualError_/objectiveNorm_;
     2891  errorCheck=maximumDualError/objectiveNorm_;
    25662892  if ((1.0-move)*errorCheck>dualTolerance()) {
    25672893    if (nextGap<gammad*(1.0-move)*errorCheck) {
     
    26462972  double largestDiagonal=0.0;
    26472973  double smallestDiagonal=1.0e50;
    2648   double * zVec = zVec_;
    2649   double * tVec = tVec_;
    2650   double * primal = solution_;
    2651   double * dual = dj_;
    2652   double * lower = lower_;
    2653   double * upper = upper_;
    2654   double * lowerSlack = lowerSlack_;
    2655   double * upperSlack = upperSlack_;
    2656   double * deltaZ = deltaZ_;
    2657   double * deltaT = deltaT_;
    2658   double * diagonal = diagonal_;
    2659   //direction vector in deltaX
    2660   double * deltaX = deltaX_;
    2661   double * cost = cost_;
    26622974  double largeGap2 = max(1.0e7,1.0e2*solutionNorm_);
     2975  //largeGap2 = 1.0e9;
    26632976  // When to start looking at killing (factor0
    26642977  double killFactor;
    2665   if (numberIterations_<50) {
    2666     killFactor = 1.0;
    2667   } else if (numberIterations_<100) {
    2668     killFactor = 10.0;
    2669     stepLength_=0.9995;
    2670   }else if (numberIterations_<150) {
    2671     killFactor = 100.0;
    2672     stepLength_=0.99995;
     2978#ifndef NO_RTTI
     2979  ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
     2980#else
     2981  ClpQuadraticObjective * quadraticObj = NULL;
     2982  if (objective_->type()==2)
     2983    quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
     2984#endif
     2985  if (!quadraticObj||1) {
     2986    if (numberIterations_<50) {
     2987      killFactor = 1.0;
     2988    } else if (numberIterations_<100) {
     2989      killFactor = 10.0;
     2990      stepLength_=max(stepLength_,0.9995);
     2991    }else if (numberIterations_<150) {
     2992      killFactor = 100.0;
     2993      stepLength_=max(stepLength_,0.99995);
     2994    } else {
     2995      killFactor = 1.0e5;
     2996      stepLength_=max(stepLength_,0.999995);
     2997    }
    26732998  } else {
    2674     killFactor = 1.0e5;
    2675     stepLength_=0.999995;
     2999    killFactor=1.0;
    26763000  }
    26773001  // put next primal into deltaSL_
    26783002  int iColumn;
    26793003  for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
    2680     if (!flagged(iColumn)) {
    2681       double thisWeight=deltaX[iColumn];
    2682       double newPrimal=primal[iColumn]+1.0*actualPrimalStep_*thisWeight;
    2683       deltaSL_[iColumn]=newPrimal;
    2684     }
     3004    double thisWeight=deltaX_[iColumn];
     3005    double newPrimal=solution_[iColumn]+1.0*actualPrimalStep_*thisWeight;
     3006    deltaSL_[iColumn]=newPrimal;
    26853007  }
    26863008  // do reduced costs
    26873009  CoinMemcpyN(dual_,numberRows_,dj_+numberColumns_);
    26883010  CoinMemcpyN(cost_,numberColumns_,dj_);
     3011  double quadraticOffset=quadraticDjs(dj_,deltaSL_,1.0);
     3012  // Save modified costs for fixed variables
     3013  CoinMemcpyN(dj_,numberColumns_,deltaSU_);
    26893014  matrix_->transposeTimes(-1.0,dual_,dj_);
    2690   double quadraticOffset=quadraticDjs(dj_,deltaSL_,1.0);
    26913015  double gamma2 = gamma_*gamma_; // gamma*gamma will be added to diagonal
    26923016  double gammaOffset=0.0;
    26933017  for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
    26943018    if (!flagged(iColumn)) {
    2695       double reducedCost=dual[iColumn];
     3019      double reducedCost=dj_[iColumn];
    26963020      bool thisKilled=false;
    2697       double zValue = zVec[iColumn] + actualDualStep_*deltaZ[iColumn];
    2698       double tValue = tVec[iColumn] + actualDualStep_*deltaT[iColumn];
    2699       zVec[iColumn]=zValue;
    2700       tVec[iColumn]=tValue;
    2701       double thisWeight=deltaX[iColumn];
    2702       double oldPrimal = primal[iColumn];
    2703       double newPrimal=primal[iColumn]+actualPrimalStep_*thisWeight;
     3021      double zValue = zVec_[iColumn] + actualDualStep_*deltaZ_[iColumn];
     3022      double wValue = wVec_[iColumn] + actualDualStep_*deltaW_[iColumn];
     3023      zVec_[iColumn]=zValue;
     3024      wVec_[iColumn]=wValue;
     3025      double thisWeight=deltaX_[iColumn];
     3026      double oldPrimal = solution_[iColumn];
     3027      double newPrimal=solution_[iColumn]+actualPrimalStep_*thisWeight;
    27043028      double dualObjectiveThis=0.0;
    27053029      double sUpper=extra;
     
    27333057        if (fakeNewBounds) {
    27343058          lower_[iColumn]=newPrimal-largeGap2;
    2735           lowerSlack[iColumn] = largeGap2;
     3059          lowerSlack_[iColumn] = largeGap2;
    27363060          upper_[iColumn]=newPrimal+largeGap2;
    2737           upperSlack[iColumn] = largeGap2;
     3061          upperSlack_[iColumn] = largeGap2;
    27383062        } else {
    27393063          lower_[iColumn]=trueLower;
    27403064          setLowerBound(iColumn);
    2741           lowerSlack[iColumn] = max(newPrimal-trueLower,1.0);
     3065          lowerSlack_[iColumn] = max(newPrimal-trueLower,1.0);
    27423066          upper_[iColumn]=trueUpper;
    27433067          setUpperBound(iColumn);
    2744           upperSlack[iColumn] = max(trueUpper-newPrimal,1.0);
     3068          upperSlack_[iColumn] = max(trueUpper-newPrimal,1.0);
    27453069        }
    27463070      } else if (fakeNewBounds) {
    27473071        lower_[iColumn]=newPrimal-largeGap2;
    2748         lowerSlack[iColumn] = largeGap2;
     3072        lowerSlack_[iColumn] = largeGap2;
    27493073        upper_[iColumn]=newPrimal+largeGap2;
    2750         upperSlack[iColumn] = largeGap2;
     3074        upperSlack_[iColumn] = largeGap2;
    27513075        // so we can just have one test
    27523076        fakeOldBounds=true;
     
    27553079      double upperBoundInfeasibility=0.0;
    27563080      if (lowerBound(iColumn)) {
    2757         double oldSlack = lowerSlack[iColumn];
     3081        double oldSlack = lowerSlack_[iColumn];
    27583082        double newSlack;
    27593083        newSlack=
    2760           lowerSlack[iColumn]+actualPrimalStep_*(oldPrimal-oldSlack
    2761                                                  + thisWeight-lower[iColumn]);
     3084          lowerSlack_[iColumn]+actualPrimalStep_*(oldPrimal-oldSlack
     3085                                                 + thisWeight-lower_[iColumn]);
    27623086        if (fakeOldBounds)
    2763           newSlack = lowerSlack[iColumn];
     3087          newSlack = lowerSlack_[iColumn];
    27643088        double epsilon = fabs(newSlack)*epsilonBase;
    27653089        if (epsilon>1.0e-5) {
     
    27723096          zValue=epsilon;
    27733097        }
    2774         double feasibleSlack=newPrimal-lower[iColumn];
     3098        double feasibleSlack=newPrimal-lower_[iColumn];
    27753099        if (feasibleSlack>0.0&&newSlack>0.0) {
    27763100          double smallGap2=smallGap;
     
    27883112          }
    27893113        }
    2790         if (zVec[iColumn]>dualTolerance) {
    2791           dualObjectiveThis+=lower[iColumn]*zVec[iColumn];
    2792         }
    2793         lowerSlack[iColumn]=newSlack;
     3114        if (zVec_[iColumn]>dualTolerance) {
     3115          dualObjectiveThis+=lower_[iColumn]*zVec_[iColumn];
     3116        }
     3117        lowerSlack_[iColumn]=newSlack;
    27943118        if (newSlack<smallerSlack) {
    27953119          smallerSlack=newSlack;
    27963120        }
    2797         lowerBoundInfeasibility = fabs(newPrimal-lowerSlack[iColumn]-lower[iColumn]);
    2798         if (lowerSlack[iColumn]<=kill*killFactor&&fabs(newPrimal-lower[iColumn])<=kill*killFactor) {
     3121        lowerBoundInfeasibility = fabs(newPrimal-lowerSlack_[iColumn]-lower_[iColumn]);
     3122        if (lowerSlack_[iColumn]<=kill*killFactor&&fabs(newPrimal-lower_[iColumn])<=kill*killFactor) {
    27993123          double step = min(actualPrimalStep_*1.1,1.0);
    2800           double newPrimal2=primal[iColumn]+step*thisWeight;
     3124          double newPrimal2=solution_[iColumn]+step*thisWeight;
    28013125          if (newPrimal2<newPrimal&&dj_[iColumn]>1.0e-5&&numberIterations_>50-40) {
    2802             newPrimal=lower[iColumn];
    2803             lowerSlack[iColumn]=0.0;
     3126            newPrimal=lower_[iColumn];
     3127            lowerSlack_[iColumn]=0.0;
    28043128            //printf("fixing %d to lower\n",iColumn);
    28053129          }
    28063130        }
    2807         if (lowerSlack[iColumn]<=kill&&fabs(newPrimal-lower[iColumn])<=kill) {
     3131        if (lowerSlack_[iColumn]<=kill&&fabs(newPrimal-lower_[iColumn])<=kill) {
    28083132          //may be better to leave at value?
    2809           newPrimal=lower[iColumn];
    2810           lowerSlack[iColumn]=0.0;
     3133          newPrimal=lower_[iColumn];
     3134          lowerSlack_[iColumn]=0.0;
    28113135          thisKilled=true;
    2812           //cout<<j<<" l "<<reducedCost<<" "<<zVec[iColumn]<<endl;
     3136          //cout<<j<<" l "<<reducedCost<<" "<<zVec_[iColumn]<<endl;
    28133137        } else {
    2814           sLower+=lowerSlack[iColumn];
     3138          sLower+=lowerSlack_[iColumn];
    28153139        }
    28163140      }
    28173141      if (upperBound(iColumn)) {
    2818         double oldSlack = upperSlack[iColumn];
     3142        double oldSlack = upperSlack_[iColumn];
    28193143        double newSlack;
    28203144        newSlack=
    2821           upperSlack[iColumn]+actualPrimalStep_*(-oldPrimal-oldSlack
    2822                                                  - thisWeight+upper[iColumn]);
     3145          upperSlack_[iColumn]+actualPrimalStep_*(-oldPrimal-oldSlack
     3146                                                 - thisWeight+upper_[iColumn]);
    28233147        if (fakeOldBounds)
    2824           newSlack = upperSlack[iColumn];
     3148          newSlack = upperSlack_[iColumn];
    28253149        double epsilon = fabs(newSlack)*epsilonBase;
    28263150        if (epsilon>1.0e-5) {
     
    28303154        //make sure reasonable
    28313155        //epsilon=1.0e-14;
    2832         if (tValue<epsilon) {
    2833           tValue=epsilon;
    2834         }
    2835         double feasibleSlack=upper[iColumn]-newPrimal;
     3156        if (wValue<epsilon) {
     3157          wValue=epsilon;
     3158        }
     3159        double feasibleSlack=upper_[iColumn]-newPrimal;
    28363160        if (feasibleSlack>0.0&&newSlack>0.0) {
    28373161          double smallGap2=smallGap;
     
    28493173          }
    28503174        }
    2851         if (tVec[iColumn]>dualTolerance) {
    2852           dualObjectiveThis-=upper[iColumn]*tVec[iColumn];
    2853         }
    2854         upperSlack[iColumn]=newSlack;
     3175        if (wVec_[iColumn]>dualTolerance) {
     3176          dualObjectiveThis-=upper_[iColumn]*wVec_[iColumn];
     3177        }
     3178        upperSlack_[iColumn]=newSlack;
    28553179        if (newSlack<smallerSlack) {
    28563180          smallerSlack=newSlack;
    28573181        }
    2858         upperBoundInfeasibility = fabs(newPrimal+upperSlack[iColumn]-upper[iColumn]);
    2859         if (upperSlack[iColumn]<=kill*killFactor&&fabs(newPrimal-upper[iColumn])<=kill*killFactor) {
     3182        upperBoundInfeasibility = fabs(newPrimal+upperSlack_[iColumn]-upper_[iColumn]);
     3183        if (upperSlack_[iColumn]<=kill*killFactor&&fabs(newPrimal-upper_[iColumn])<=kill*killFactor) {
    28603184          double step = min(actualPrimalStep_*1.1,1.0);
    2861           double newPrimal2=primal[iColumn]+step*thisWeight;
     3185          double newPrimal2=solution_[iColumn]+step*thisWeight;
    28623186          if (newPrimal2>newPrimal&&dj_[iColumn]<-1.0e-5&&numberIterations_>50-40) {
    2863             newPrimal=upper[iColumn];
    2864             upperSlack[iColumn]=0.0;
     3187            newPrimal=upper_[iColumn];
     3188            upperSlack_[iColumn]=0.0;
    28653189            //printf("fixing %d to upper\n",iColumn);
    28663190          }
    28673191        }
    2868         if (upperSlack[iColumn]<=kill&&fabs(newPrimal-upper[iColumn])<=kill) {
     3192        if (upperSlack_[iColumn]<=kill&&fabs(newPrimal-upper_[iColumn])<=kill) {
    28693193          //may be better to leave at value?
    2870           newPrimal=upper[iColumn];
    2871           upperSlack[iColumn]=0.0;
     3194          newPrimal=upper_[iColumn];
     3195          upperSlack_[iColumn]=0.0;
    28723196          thisKilled=true;
    28733197        } else {
    2874           sUpper+=upperSlack[iColumn];
    2875         }
    2876       }
    2877       primal[iColumn]=newPrimal;
     3198          sUpper+=upperSlack_[iColumn];
     3199        }
     3200      }
     3201      solution_[iColumn]=newPrimal;
    28783202      if (fabs(newPrimal)>solutionNorm) {
    28793203        solutionNorm=fabs(newPrimal);
     
    28863210        }
    28873211        double dualInfeasibility=
    2888           reducedCost-zVec[iColumn]+tVec[iColumn]+gammaTerm*newPrimal;
     3212          reducedCost-zVec_[iColumn]+wVec_[iColumn]+gammaTerm*newPrimal;
    28893213        if (fabs(dualInfeasibility)>dualTolerance) {
    28903214#if 0
     
    28923216            // To improve we could reduce t and/or increase z
    28933217            if (lowerBound(iColumn)) {
    2894               double complementarity =zVec_[iColumn]*lowerSlack[iColumn];
     3218              double complementarity =zVec_[iColumn]*lowerSlack_[iColumn];
    28953219              if (complementarity<nextMu) {
    28963220                double change=
    28973221                  min(dualInfeasibility,
    2898                       (nextMu-complementarity)/lowerSlack[iColumn]);
     3222                      (nextMu-complementarity)/lowerSlack_[iColumn]);
    28993223                dualInfeasibility -= change;
    29003224                printf("%d lb locomp %g - dual inf from %g to %g\n",
     
    29063230            }
    29073231            if (upperBound(iColumn)) {
    2908               double complementarity =tVec_[iColumn]*upperSlack[iColumn];
     3232              double complementarity =wVec_[iColumn]*upperSlack_[iColumn];
    29093233              if (complementarity>nextMu) {
    29103234                double change=
    29113235                  min(dualInfeasibility,
    2912                       (complementarity-nextMu)/upperSlack[iColumn]);
     3236                      (complementarity-nextMu)/upperSlack_[iColumn]);
    29133237                dualInfeasibility -= change;
    29143238                printf("%d ub hicomp %g - dual inf from %g to %g\n",
    29153239                       iColumn,complementarity,dualInfeasibility+change,
    29163240                       dualInfeasibility);
    2917                 tVec_[iColumn] -= change;
    2918                 tValue = max(tVec_[iColumn],1.0e-12);
     3241                wVec_[iColumn] -= change;
     3242                wValue = max(wVec_[iColumn],1.0e-12);
    29193243              }
    29203244            }
     
    29223246            // To improve we could reduce z and/or increase t
    29233247            if (lowerBound(iColumn)) {
    2924               double complementarity =zVec_[iColumn]*lowerSlack[iColumn];
     3248              double complementarity =zVec_[iColumn]*lowerSlack_[iColumn];
    29253249              if (complementarity>nextMu) {
    29263250                double change=
    29273251                  max(dualInfeasibility,
    2928                       (nextMu-complementarity)/lowerSlack[iColumn]);
     3252                      (nextMu-complementarity)/lowerSlack_[iColumn]);
    29293253                dualInfeasibility -= change;
    29303254                printf("%d lb hicomp %g - dual inf from %g to %g\n",
     
    29363260            }
    29373261            if (upperBound(iColumn)) {
    2938               double complementarity =tVec_[iColumn]*upperSlack[iColumn];
     3262              double complementarity =wVec_[iColumn]*upperSlack_[iColumn];
    29393263              if (complementarity<nextMu) {
    29403264                double change=
    29413265                  max(dualInfeasibility,
    2942                       (complementarity-nextMu)/upperSlack[iColumn]);
     3266                      (complementarity-nextMu)/upperSlack_[iColumn]);
    29433267                dualInfeasibility -= change;
    29443268                printf("%d ub locomp %g - dual inf from %g to %g\n",
    29453269                       iColumn,complementarity,dualInfeasibility+change,
    29463270                       dualInfeasibility);
    2947                 tVec_[iColumn] -= change;
    2948                 tValue = max(tVec_[iColumn],1.0e-12);
     3271                wVec_[iColumn] -= change;
     3272                wValue = max(wVec_[iColumn],1.0e-12);
    29493273              }
    29503274            }
     
    29623286        if (dualInfeasibility>maximumDualError) {
    29633287          //printf("bad dual %d %g\n",iColumn,
    2964           // reducedCost-zVec[iColumn]+tVec[iColumn]+gammaTerm*newPrimal);
     3288          // reducedCost-zVec_[iColumn]+wVec_[iColumn]+gammaTerm*newPrimal);
    29653289          maximumDualError=dualInfeasibility;
    29663290        }
     
    29743298        }
    29753299#if 1
    2976         double divisor = sLower*tValue+sUpper*zValue+gammaTerm*sLower*sUpper;
     3300        double divisor = sLower*wValue+sUpper*zValue+gammaTerm*sLower*sUpper;
    29773301        double diagonalValue=(sUpper*sLower)/divisor;
    29783302#else
    2979         double divisor = sLower*tValue+sUpper*zValue+gammaTerm*sLower*sUpper;
     3303        double divisor = sLower*wValue+sUpper*zValue+gammaTerm*sLower*sUpper;
    29803304        double diagonalValue2=(sUpper*sLower)/divisor;
    29813305        double diagonalValue;
    29823306        if (!lowerBound(iColumn)) {
    2983           diagonalValue = tValue/sUpper + gammaTerm;
     3307          diagonalValue = wValue/sUpper + gammaTerm;
    29843308        } else if (!upperBound(iColumn)) {
    29853309          diagonalValue = zValue/sLower + gammaTerm;
    29863310        } else {
    2987           diagonalValue = zValue/sLower + tValue/sUpper + gammaTerm;
     3311          diagonalValue = zValue/sLower + wValue/sUpper + gammaTerm;
    29883312        }
    29893313        diagonalValue = 1.0/diagonalValue;
    29903314#endif
    2991         diagonal[iColumn]=diagonalValue;
     3315        diagonal_[iColumn]=diagonalValue;
    29923316        //FUDGE
    29933317        if (diagonalValue>diagonalLimit) {
    29943318          std::cout<<"large diagonal "<<diagonalValue<<std::endl;
    2995           diagonal[iColumn]=diagonalLimit;
     3319          diagonal_[iColumn]=diagonalLimit;
    29963320        }
    29973321        if (diagonalValue<1.0e-10) {
     
    30043328          smallestDiagonal=diagonalValue;
    30053329        }
    3006         deltaX[iColumn]=0.0;
     3330        deltaX_[iColumn]=0.0;
    30073331      } else {
    30083332        numberKilled++;
    3009         diagonal[iColumn]=0.0;
    3010         zVec[iColumn]=0.0;
    3011         tVec[iColumn]=0.0;
     3333        diagonal_[iColumn]=0.0;
     3334        zVec_[iColumn]=0.0;
     3335        wVec_[iColumn]=0.0;
    30123336        setFlagged(iColumn);
    30133337        setFixedOrFree(iColumn);
    3014         deltaX[iColumn]=newPrimal;
    3015         offsetObjective+=newPrimal*cost[iColumn];
     3338        deltaX_[iColumn]=newPrimal;
     3339        offsetObjective+=newPrimal*deltaSU_[iColumn];
    30163340      }
    30173341    } else {
    3018       deltaX[iColumn]=primal[iColumn];
    3019       diagonal[iColumn]=0.0;
    3020       offsetObjective+=primal[iColumn]*cost[iColumn];
    3021       if (upper[iColumn]-lower[iColumn]>1.0e-5) {
    3022         if (primal[iColumn]<lower[iColumn]+1.0e-8&&dual[iColumn]<-1.0e-8) {
    3023           if (-dual[iColumn]>maximumDJInfeasibility) {
    3024             maximumDJInfeasibility=-dual[iColumn];
     3342      deltaX_[iColumn]=solution_[iColumn];
     3343      diagonal_[iColumn]=0.0;
     3344      offsetObjective+=solution_[iColumn]*deltaSU_[iColumn];
     3345      if (upper_[iColumn]-lower_[iColumn]>1.0e-5) {
     3346        if (solution_[iColumn]<lower_[iColumn]+1.0e-8&&dj_[iColumn]<-1.0e-8) {
     3347          if (-dj_[iColumn]>maximumDJInfeasibility) {
     3348            maximumDJInfeasibility=-dj_[iColumn];
    30253349          }
    30263350        }
    3027         if (primal[iColumn]>upper[iColumn]-1.0e-8&&dual[iColumn]>1.0e-8) {
    3028           if (dual[iColumn]>maximumDJInfeasibility) {
    3029             maximumDJInfeasibility=dual[iColumn];
     3351        if (solution_[iColumn]>upper_[iColumn]-1.0e-8&&dj_[iColumn]>1.0e-8) {
     3352          if (dj_[iColumn]>maximumDJInfeasibility) {
     3353            maximumDJInfeasibility=dj_[iColumn];
    30303354          }
    30313355        }
    30323356      }
    30333357    }
    3034     primalObjectiveValue+=primal[iColumn]*cost[iColumn];
     3358    primalObjectiveValue+=solution_[iColumn]*cost_[iColumn];
    30353359  }
    30363360  handler_->message(CLP_BARRIER_DIAGONAL,messages_)
     
    30483372#endif
    30493373  // update rhs region
    3050   multiplyAdd(deltaX+numberColumns_,numberRows_,-1.0,rhsFixRegion_,1.0);
    3051   matrix_->times(1.0,deltaX,rhsFixRegion_);
     3374  multiplyAdd(deltaX_+numberColumns_,numberRows_,-1.0,rhsFixRegion_,1.0);
     3375  matrix_->times(1.0,deltaX_,rhsFixRegion_);
    30523376  primalObjectiveValue += 0.5*gamma2*gammaOffset+0.5*quadraticOffset;
    30533377  if (quadraticOffset) {
    3054     //printf("gamma offset %g %g, quadoffset %g\n",gammaOffset,gamma2*gammaOffset,quadraticOffset);
     3378    //  printf("gamma offset %g %g, quadoffset %g\n",gammaOffset,gamma2*gammaOffset,quadraticOffset);
    30553379  }
    30563380 
     
    31873511  return numberKilled;
    31883512}
    3189 //  Save info on products of affine deltaT*deltaW and deltaS*deltaZ
     3513//  Save info on products of affine deltaW*deltaW and deltaS*deltaZ
    31903514double
    31913515ClpPredictorCorrector::affineProduct()
    31923516{
    3193   double * deltaZ = deltaZ_;
    3194   double * deltaT = deltaT_;
    3195   double * lower = lower_;
    3196   double * upper = upper_;
    3197   double * lowerSlack = lowerSlack_;
    3198   double * upperSlack = upperSlack_;
    3199   double * primal = solution_;
    3200   //direction vector in deltaX
    3201   double * deltaX = deltaX_;
    32023517  double product = 0.0;
    32033518  //IF zVec starts as 0 then deltaZ always zero
     
    32063521  //out all we want
    32073522  for (int iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
    3208     double w3=deltaZ[iColumn]*deltaX[iColumn];
    3209     double w4=-deltaT[iColumn]*deltaX[iColumn];
     3523    double w3=deltaZ_[iColumn]*deltaX_[iColumn];
     3524    double w4=-deltaW_[iColumn]*deltaX_[iColumn];
    32103525    if (lowerBound(iColumn)) {
    3211       w3+=deltaZ[iColumn]*(primal[iColumn]-lowerSlack[iColumn]-lower[iColumn]);
     3526      w3+=deltaZ_[iColumn]*(solution_[iColumn]-lowerSlack_[iColumn]-lower_[iColumn]);
    32123527      product+=w3;
    32133528    }
    32143529    if (upperBound(iColumn)) {
    3215       w4+=deltaT[iColumn]*(-primal[iColumn]-upperSlack[iColumn]+upper[iColumn]);
     3530      w4+=deltaW_[iColumn]*(-solution_[iColumn]-upperSlack_[iColumn]+upper_[iColumn]);
    32163531      product+=w4;
    32173532    }
     
    32233538ClpPredictorCorrector::debugMove(int phase,double primalStep, double dualStep)
    32243539{
     3540#ifndef SOME_DEBUG
    32253541  return;
    3226   double * zVec = zVec_;
    3227   double * tVec = tVec_;
    3228   double * primal = solution_;
    3229   double * lower = lower_;
    3230   double * upper = upper_;
    3231   double * lowerSlack = lowerSlack_;
    3232   double * upperSlack = upperSlack_;
    3233   double * deltaZ = deltaZ_;
    3234   double * deltaT = deltaT_;
    3235   //direction vector in deltaX
    3236   double * deltaX = deltaX_;
    3237   double * cost = cost_;
     3542#endif
    32383543  int numberTotal = numberRows_+numberColumns_;
    32393544  double * dualNew = ClpCopyOfArray(dual_,numberRows_);
     
    32483553  CoinMemcpyN(cost_,numberColumns_,djNew);
    32493554  matrix_->transposeTimes(-1.0,dualNew,djNew);
    3250   double quadraticOffset=quadraticDjs(djNew,primalNew,1.0);
    32513555  // update x
    32523556  int iColumn;
    32533557  for (iColumn=0;iColumn<numberTotal;iColumn++) {
    32543558    if (!flagged(iColumn))
    3255       primalNew[iColumn] +=primalStep*deltaX[iColumn];
    3256   }
     3559      primalNew[iColumn] +=primalStep*deltaX_[iColumn];
     3560  }
     3561  double quadraticOffset=quadraticDjs(djNew,primalNew,1.0);
    32573562  CoinZeroN(errorRegionNew,numberRows_);
    32583563  CoinZeroN(rhsFixRegionNew,numberRows_);
     
    32813586    if (!flagged(iColumn)) {
    32823587      double reducedCost=djNew[iColumn];
    3283       double zValue = zVec[iColumn] + dualStep*deltaZ[iColumn];
    3284       double tValue = tVec[iColumn] + dualStep*deltaT[iColumn];
    3285       double thisWeight=deltaX[iColumn];
     3588      double zValue = zVec_[iColumn] + dualStep*deltaZ_[iColumn];
     3589      double wValue = wVec_[iColumn] + dualStep*deltaW_[iColumn];
     3590      double thisWeight=deltaX_[iColumn];
    32863591      double oldPrimal = solution_[iColumn];
    32873592      double newPrimal=primalNew[iColumn];
     
    32893594      double upperBoundInfeasibility=0.0;
    32903595      if (lowerBound(iColumn)) {
    3291         double oldSlack = lowerSlack[iColumn];
     3596        double oldSlack = lowerSlack_[iColumn];
    32923597        double newSlack=
    3293           lowerSlack[iColumn]+primalStep*(oldPrimal-oldSlack
    3294                                                  + thisWeight-lower[iColumn]);
     3598          lowerSlack_[iColumn]+primalStep*(oldPrimal-oldSlack
     3599                                                 + thisWeight-lower_[iColumn]);
    32953600        if (zValue>dualTolerance) {
    3296           dualObjectiveValue+=lower[iColumn]*zVec[iColumn];
    3297         }
    3298         lowerBoundInfeasibility = fabs(newPrimal-newSlack-lower[iColumn]);
     3601          dualObjectiveValue+=lower_[iColumn]*zVec_[iColumn];
     3602        }
     3603        lowerBoundInfeasibility = fabs(newPrimal-newSlack-lower_[iColumn]);
    32993604        newGap += newSlack*zValue;
    33003605      }
    33013606      if (upperBound(iColumn)) {
    3302         double oldSlack = upperSlack[iColumn];
     3607        double oldSlack = upperSlack_[iColumn];
    33033608        double newSlack=
    3304           upperSlack[iColumn]+primalStep*(-oldPrimal-oldSlack
    3305                                                  - thisWeight+upper[iColumn]);
    3306         if (tValue>dualTolerance) {
    3307           dualObjectiveValue-=upper[iColumn]*tVec[iColumn];
    3308         }
    3309         upperBoundInfeasibility = fabs(newPrimal+newSlack-upper[iColumn]);
    3310         newGap += newSlack*tValue;
     3609          upperSlack_[iColumn]+primalStep*(-oldPrimal-oldSlack
     3610                                                 - thisWeight+upper_[iColumn]);
     3611        if (wValue>dualTolerance) {
     3612          dualObjectiveValue-=upper_[iColumn]*wVec_[iColumn];
     3613        }
     3614        upperBoundInfeasibility = fabs(newPrimal+newSlack-upper_[iColumn]);
     3615        newGap += newSlack*wValue;
    33113616      }
    33123617      if (fabs(newPrimal)>solutionNorm) {
     
    33193624      }
    33203625      double dualInfeasibility=
    3321         reducedCost-zValue+tValue+gammaTerm*newPrimal;
     3626        reducedCost-zValue+wValue+gammaTerm*newPrimal;
    33223627      if (fabs(dualInfeasibility)>dualTolerance) {
    33233628        dualFake+=newPrimal*dualInfeasibility;
     
    33323637      if (dualInfeasibility>maximumDualError) {
    33333638        //printf("bad dual %d %g\n",iColumn,
    3334         // reducedCost-zVec[iColumn]+tVec[iColumn]+gammaTerm*newPrimal);
     3639        // reducedCost-zVec_[iColumn]+wVec_[iColumn]+gammaTerm*newPrimal);
    33353640        maximumDualError=dualInfeasibility;
    33363641      }
     
    33383643      djNew[iColumn]=0.0;
    33393644    } else {
    3340       offsetObjective+=primalNew[iColumn]*cost[iColumn];
    3341       if (upper[iColumn]-lower[iColumn]>1.0e-5) {
    3342         if (primalNew[iColumn]<lower[iColumn]+1.0e-8&&djNew[iColumn]<-1.0e-8) {
     3645      offsetObjective+=primalNew[iColumn]*cost_[iColumn];
     3646      if (upper_[iColumn]-lower_[iColumn]>1.0e-5) {
     3647        if (primalNew[iColumn]<lower_[iColumn]+1.0e-8&&djNew[iColumn]<-1.0e-8) {
    33433648          if (-djNew[iColumn]>maximumDjInfeasibility) {
    33443649            maximumDjInfeasibility=-djNew[iColumn];
    33453650          }
    33463651        }
    3347         if (primalNew[iColumn]>upper[iColumn]-1.0e-8&&djNew[iColumn]>1.0e-8) {
     3652        if (primalNew[iColumn]>upper_[iColumn]-1.0e-8&&djNew[iColumn]>1.0e-8) {
    33483653          if (djNew[iColumn]>maximumDjInfeasibility) {
    33493654            maximumDjInfeasibility=djNew[iColumn];
     
    33533658      djNew[iColumn]=primalNew[iColumn];
    33543659    }
    3355     primalObjectiveValue+=primal[iColumn]*cost[iColumn];
     3660    primalObjectiveValue+=solution_[iColumn]*cost_[iColumn];
    33563661  }
    33573662  // update rhs region
  • trunk/ClpPresolve.cpp

    r378 r393  
    1313#include "CoinPackedMatrix.hpp"
    1414#include "ClpSimplex.hpp"
     15#include "ClpQuadraticObjective.hpp"
    1516
    1617#include "ClpPresolve.hpp"
     
    215216  }
    216217  // Now check solution
     218  double offset;
    217219  memcpy(originalModel_->dualColumnSolution(),
    218          originalModel_->objective(),ncols_*sizeof(double));
     220         originalModel_->objectiveAsObject()->gradient(originalModel_,
     221                                                       originalModel_->primalColumnSolution(),
     222                                                       offset,true),ncols_*sizeof(double));
    219223  originalModel_->transposeTimes(-1.0,
    220224                                 originalModel_->dualRowSolution(),
     
    371375  // later just do individually
    372376  bool doDualStuff = (presolvedModel_->integerInformation()==NULL);
     377  if (prob->anyProhibited())
     378    doDualStuff=false;
    373379
    374380#if     CHECK_CONSISTENCY
     
    862868  ClpDisjointCopyN(si->getColLower(), ncols, clo_);
    863869  ClpDisjointCopyN(si->getColUpper(), ncols, cup_);
    864   ClpDisjointCopyN(si->getObjCoefficients(), ncols, cost_);
     870  //ClpDisjointCopyN(si->getObjCoefficients(), ncols, cost_);
     871  double offset;
     872  ClpDisjointCopyN(si->objectiveAsObject()->gradient(si,si->getColSolution(),offset,true), ncols, cost_);
    865873  ClpDisjointCopyN(si->getRowLower(), nrows,  rlo_);
    866874  ClpDisjointCopyN(si->getRowUpper(), nrows,  rup_);
     
    991999  }
    9921000
     1001#ifndef NO_RTTI
     1002    ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(si->objectiveAsObject()));
     1003#else
     1004    ClpQuadraticObjective * quadraticObj = NULL;
     1005    if (si->objectiveAsObject()->type()==2)
     1006      quadraticObj = (static_cast< ClpQuadraticObjective*>(si->objectiveAsObject()));
     1007#endif
    9931008  // Set up prohibited bits if needed
    9941009  if (nonLinearValue) {
     
    10081023        setColProhibited(icol);
    10091024    }
     1025  } else if (quadraticObj) {
     1026    CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
     1027    //const int * columnQuadratic = quadratic->getIndices();
     1028    //const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
     1029    const int * columnQuadraticLength = quadratic->getVectorLengths();
     1030    //double * quadraticElement = quadratic->getMutableElements();
     1031    int numberColumns = quadratic->getNumCols();
     1032    anyProhibited_ = true;
     1033    for (int iColumn=0;iColumn<numberColumns;iColumn++) {
     1034      if (columnQuadraticLength[iColumn]) {
     1035        setColProhibited(iColumn);
     1036        //printf("%d prohib\n",iColumn);
     1037      }
     1038    }
    10101039  } else {
    10111040    anyProhibited_ = false;
     
    10621091  si->loadProblem(ncols_, nrows_, mcstrt_, hrow_, colels_, hincol_,
    10631092                 clo_, cup_, cost_, rlo_, rup_);
    1064 
    10651093  delete [] si->integerInformation();
    10661094  int numberIntegers=0;
     
    12911319
    12921320    // Do presolve
    1293 
    12941321    paction_ = presolve(&prob);
    12951322
     
    13511378      int ncolsNow = presolvedModel_->getNumCols();
    13521379      memcpy(originalColumn_,prob.originalColumn_,ncolsNow*sizeof(int));
     1380#ifndef NO_RTTI
     1381      ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(originalModel->objectiveAsObject()));
     1382#else
     1383      ClpQuadraticObjective * quadraticObj = NULL;
     1384      if (originalModel->objectiveAsObject()->type()==2)
     1385        quadraticObj = (static_cast< ClpQuadraticObjective*>(originalModel->objectiveAsObject()));
     1386#endif
     1387      if (quadraticObj) {
     1388        // set up for subset
     1389        char * mark = new char [ncols_];
     1390        memset(mark,0,ncols_);
     1391        CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
     1392        //const int * columnQuadratic = quadratic->getIndices();
     1393        //const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
     1394        const int * columnQuadraticLength = quadratic->getVectorLengths();
     1395        //double * quadraticElement = quadratic->getMutableElements();
     1396        int numberColumns = quadratic->getNumCols();
     1397        ClpQuadraticObjective * newObj = new ClpQuadraticObjective(*quadraticObj,
     1398                                                                   ncolsNow,
     1399                                                                   originalColumn_);
     1400        // and modify linear and check
     1401        double * linear = newObj->linearObjective();
     1402        memcpy(linear,presolvedModel_->objective(),ncolsNow*sizeof(double));
     1403        int iColumn;
     1404        for ( iColumn=0;iColumn<numberColumns;iColumn++) {
     1405          if (columnQuadraticLength[iColumn])
     1406            mark[iColumn]=1;
     1407        }
     1408        // and new
     1409        quadratic = newObj->quadraticObjective();
     1410        columnQuadraticLength = quadratic->getVectorLengths();
     1411        int numberColumns2 = quadratic->getNumCols();
     1412        for ( iColumn=0;iColumn<numberColumns2;iColumn++) {
     1413          if (columnQuadraticLength[iColumn])
     1414            mark[originalColumn_[iColumn]]=0;
     1415        }
     1416        presolvedModel_->setObjective(newObj);
     1417        delete newObj;
     1418        // final check
     1419        for ( iColumn=0;iColumn<numberColumns;iColumn++)
     1420          if (mark[iColumn])
     1421            printf("Quadratic column %d modified - may be okay\n",iColumn);
     1422        delete [] mark;
     1423      }
    13531424      delete [] prob.originalColumn_;
    13541425      prob.originalColumn_=NULL;
  • trunk/ClpPrimalColumnPivot.cpp

    r225 r393  
    1616ClpPrimalColumnPivot::ClpPrimalColumnPivot () :
    1717  model_(NULL),
    18   type_(-1)
     18  type_(-1),
     19  looksOptimal_(false)
    1920{
    2021
     
    2627ClpPrimalColumnPivot::ClpPrimalColumnPivot (const ClpPrimalColumnPivot & source) :
    2728  model_(source.model_),
    28   type_(source.type_)
     29  type_(source.type_),
     30  looksOptimal_(source.looksOptimal_)
    2931
    3032
     
    4850    type_ = rhs.type_;
    4951    model_ = rhs.model_;
     52    looksOptimal_ = rhs.looksOptimal_;
    5053  }
    5154  return *this;
  • trunk/ClpPrimalColumnSteepest.cpp

    r389 r393  
    33343334ClpPrimalColumnSteepest::looksOptimal() const
    33353335{
     3336  if (looksOptimal_)
     3337    return true; // user overrode
    33363338  //**** THIS MUST MATCH the action coding above
    33373339  double tolerance=model_->currentDualTolerance();
  • trunk/ClpQuadraticObjective.cpp

    r389 r393  
    44#include "CoinPragma.hpp"
    55#include "CoinHelperFunctions.hpp"
    6 #include "ClpModel.hpp"
     6#include "CoinIndexedVector.hpp"
     7#include "ClpFactorization.hpp"
     8#include "ClpSimplex.hpp"
    79#include "ClpQuadraticObjective.hpp"
    810//#############################################################################
     
    2224  numberColumns_=0;
    2325  numberExtendedColumns_=0;
     26  activated_=0;
    2427}
    2528
     
    5457  quadraticObjective_=NULL;
    5558  gradient_ = NULL;
     59  activated_=1;
    5660}
    5761
     
    192196  int extra = rhs.numberExtendedColumns_-rhs.numberColumns_;
    193197  numberColumns_=0;
     198  numberExtendedColumns_ = numberColumns+extra;
    194199  if (numberColumns>0) {
    195200    // check valid lists
     
    203208                      "ClpQuadraticObjective");
    204209    numberColumns_ = numberColumns;
    205     numberExtendedColumns_ = numberColumns+extra;
    206210    objective_ = new double[numberExtendedColumns_];
    207211    for (i=0;i<numberColumns_;i++)
     
    215219      memcpy(gradient_+numberColumns_,rhs.gradient_+rhs.numberColumns_,
    216220             (numberExtendedColumns_-numberColumns_)*sizeof(double));
    217     }
     221    } else {
     222      gradient_=NULL;
     223    }
     224  } else {
     225    gradient_=NULL;
     226    objective_=NULL;
    218227  }
    219228  if (rhs.quadraticObjective_) {
     
    272281// Returns gradient
    273282double * 
    274 ClpQuadraticObjective::gradient(const double * solution, double & offset,bool refresh)
     283ClpQuadraticObjective::gradient(const ClpSimplex * model,
     284                                const double * solution, double & offset,bool refresh,
     285                                int includeLinear)
    275286{
    276287  offset=0.0;
    277   if (!quadraticObjective_||!solution) {
    278     return objective_;
    279   } else {
     288  bool scaling=false;
     289  if (model&&(model->rowScale()||
     290              model->objectiveScale()!=1.0||model->optimizationDirection()!=1.0))
     291    scaling=true;
     292  const double * cost = NULL;
     293  if (model)
     294    cost = model->costRegion();
     295  if (!cost) {
     296    // not in solve
     297    cost = objective_;
     298    scaling=false;
     299  }
     300  if (!scaling) {
     301    if (!quadraticObjective_||!solution||!activated_) {
     302      return objective_;
     303    } else {
     304      if (refresh||!gradient_) {
     305        if (!gradient_)
     306          gradient_ = new double[numberExtendedColumns_];
     307        const int * columnQuadratic = quadraticObjective_->getIndices();
     308        const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
     309        const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
     310        const double * quadraticElement = quadraticObjective_->getElements();
     311        offset=0.0;
     312        // use current linear cost region
     313        if (includeLinear==1)
     314          memcpy(gradient_,cost,numberExtendedColumns_*sizeof(double));
     315        else if (includeLinear==2)
     316          memcpy(gradient_,objective_,numberExtendedColumns_*sizeof(double));
     317        else
     318          memset(gradient_,0,numberExtendedColumns_*sizeof(double));
     319        if (activated_) {
     320          int iColumn;
     321          for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     322            double valueI = solution[iColumn];
     323            CoinBigIndex j;
     324            for (j=columnQuadraticStart[iColumn];
     325                 j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     326              int jColumn = columnQuadratic[j];
     327              double valueJ = solution[jColumn];
     328              double elementValue = quadraticElement[j];
     329              if (iColumn!=jColumn) {
     330                offset += valueI*valueJ*elementValue;
     331                double gradientI = valueJ*elementValue;
     332                double gradientJ = valueI*elementValue;
     333                gradient_[iColumn] += gradientI;
     334                gradient_[jColumn] += gradientJ;
     335              } else {
     336                offset += 0.5*valueI*valueI*elementValue;
     337                double gradientI = valueI*elementValue;
     338                gradient_[iColumn] += gradientI;
     339              }
     340            }
     341          }
     342        }
     343      }
     344      return gradient_;
     345    }
     346  } else {
     347    // do scaling
     348    assert (solution);
    280349    if (refresh||!gradient_) {
    281350      if (!gradient_)
    282351        gradient_ = new double[numberExtendedColumns_];
     352      double direction = model->optimizationDirection()*model->objectiveScale();
     353      // direction is actually scale out not scale in
     354      if (direction)
     355        direction = 1.0/direction;
    283356      const int * columnQuadratic = quadraticObjective_->getIndices();
    284357      const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
    285358      const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
    286359      const double * quadraticElement = quadraticObjective_->getElements();
    287       double offset=0.0;
    288       memcpy(gradient_,objective_,numberExtendedColumns_*sizeof(double));
    289360      int iColumn;
    290       for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    291         double valueI = solution[iColumn];
    292         int j;
    293         for (j=columnQuadraticStart[iColumn];
    294              j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
    295           int jColumn = columnQuadratic[j];
    296           double valueJ = solution[jColumn];
    297           double elementValue = quadraticElement[j];
    298           elementValue *= 0.5;
    299           offset += valueI*valueJ*elementValue;
    300           double gradientI = valueJ*elementValue;
    301           double gradientJ = valueI*elementValue;
    302           offset -= gradientI*valueI;
    303           gradient_[iColumn] += gradientI;
    304           offset -= gradientJ*valueJ;
    305           gradient_[jColumn] += gradientJ;
     361      const double * columnScale = model->columnScale();
     362      // use current linear cost region (already scaled)
     363      if (includeLinear==1) {
     364        memcpy(gradient_,model->costRegion(),numberExtendedColumns_*sizeof(double));
     365      } else if (includeLinear==2) {
     366        memset(gradient_+numberColumns_,0,(numberExtendedColumns_-numberColumns_)*sizeof(double));
     367        if (!columnScale) {
     368          for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     369            gradient_[iColumn]= objective_[iColumn]*direction;
     370          }
     371        } else {
     372          for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     373            gradient_[iColumn]= objective_[iColumn]*direction*columnScale[iColumn];
     374          }
     375        }
     376      } else {
     377        memset(gradient_,0,numberExtendedColumns_*sizeof(double));
     378      }
     379      if (!columnScale) {
     380        if (activated_) {
     381          for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     382            double valueI = solution[iColumn];
     383            CoinBigIndex j;
     384            for (j=columnQuadraticStart[iColumn];
     385                 j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     386              int jColumn = columnQuadratic[j];
     387              double valueJ = solution[jColumn];
     388              double elementValue = quadraticElement[j];
     389              elementValue *= direction;
     390              if (iColumn!=jColumn) {
     391                offset += valueI*valueJ*elementValue;
     392                double gradientI = valueJ*elementValue;
     393                double gradientJ = valueI*elementValue;
     394                gradient_[iColumn] += gradientI;
     395                gradient_[jColumn] += gradientJ;
     396              } else {
     397                offset += 0.5*valueI*valueI*elementValue;
     398                double gradientI = valueI*elementValue;
     399                gradient_[iColumn] += gradientI;
     400              }
     401            }
     402          }
     403        }
     404      } else {
     405        // scaling
     406        if (activated_) {
     407          for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     408            double valueI = solution[iColumn];
     409            double scaleI = columnScale[iColumn]*direction;
     410            CoinBigIndex j;
     411            for (j=columnQuadraticStart[iColumn];
     412                 j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     413              int jColumn = columnQuadratic[j];
     414              double valueJ = solution[jColumn];
     415              double elementValue = quadraticElement[j];
     416              double scaleJ = columnScale[jColumn];
     417              elementValue *= scaleI*scaleJ;
     418              if (iColumn!=jColumn) {
     419                offset += valueI*valueJ*elementValue;
     420                double gradientI = valueJ*elementValue;
     421                double gradientJ = valueI*elementValue;
     422                gradient_[iColumn] += gradientI;
     423                gradient_[jColumn] += gradientJ;
     424              } else {
     425                offset += 0.5*valueI*valueI*elementValue;
     426                double gradientI = valueI*elementValue;
     427                gradient_[iColumn] += gradientI;
     428              }
     429            }
     430          }
    306431        }
    307432      }
     
    482607  quadraticObjective_ = NULL;
    483608}
     609/* Returns reduced gradient.Returns an offset (to be added to current one).
     610 */
     611double
     612ClpQuadraticObjective::reducedGradient(ClpSimplex * model, double * region,
     613                                       bool useFeasibleCosts)
     614{
     615  int numberRows = model->numberRows();
     616  int numberColumns=model->numberColumns();
     617 
     618  //work space
     619  CoinIndexedVector  * workSpace = model->rowArray(0);
     620 
     621  CoinIndexedVector arrayVector;
     622  arrayVector.reserve(numberRows+1);
     623 
     624  int iRow;
     625#ifdef CLP_DEBUG
     626  workSpace->checkClear();
     627#endif
     628  double * array = arrayVector.denseVector();
     629  int * index = arrayVector.getIndices();
     630  int number=0;
     631  const double * costNow = gradient(model,model->solutionRegion(),offset_,
     632                                    true,useFeasibleCosts ? 2 : 1);
     633  double * cost = model->costRegion();
     634  const int * pivotVariable = model->pivotVariable();
     635  for (iRow=0;iRow<numberRows;iRow++) {
     636    int iPivot=pivotVariable[iRow];
     637    double value;
     638    if (iPivot<numberColumns)
     639      value = costNow[iPivot];
     640    else if (!useFeasibleCosts)
     641      value = cost[iPivot];
     642    else
     643      value=0.0;
     644    if (value) {
     645      array[iRow]=value;
     646      index[number++]=iRow;
     647    }
     648  }
     649  arrayVector.setNumElements(number);
     650
     651  // Btran basic costs
     652  double * work = workSpace->denseVector();
     653  model->factorization()->updateColumnTranspose(workSpace,&arrayVector);
     654  ClpFillN(work,numberRows,0.0);
     655  // now look at dual solution
     656  double * rowReducedCost = region+numberColumns;
     657  double * dual = rowReducedCost;
     658  const double * rowCost = cost+numberColumns;
     659  for (iRow=0;iRow<numberRows;iRow++) {
     660    dual[iRow]=array[iRow];
     661  }
     662  double * dj = region;
     663  ClpDisjointCopyN(costNow,numberColumns,dj);
     664 
     665  model->transposeTimes(-1.0,dual,dj);
     666  for (iRow=0;iRow<numberRows;iRow++) {
     667    // slack
     668    double value = dual[iRow];
     669    value += rowCost[iRow];
     670    rowReducedCost[iRow]=value;
     671  }
     672  return offset_;
     673}
     674/* Returns step length which gives minimum of objective for
     675   solution + theta * change vector up to maximum theta.
     676   
     677   arrays are numberColumns+numberRows
     678*/
     679double
     680ClpQuadraticObjective::stepLength(ClpSimplex * model,
     681                                  const double * solution,
     682                                  const double * change,
     683                                  double maximumTheta,
     684                                  double & currentObj,
     685                                  double & predictedObj,
     686                                  double & thetaObj)
     687{
     688  const double * cost = model->costRegion();
     689  bool inSolve=true;
     690  if (!cost) {
     691    // not in solve
     692    cost = objective_;
     693    inSolve=false;
     694  }
     695  double delta=0.0;
     696  double linearCost =0.0;
     697  int numberRows = model->numberRows();
     698  int numberColumns = model->numberColumns();
     699  int numberTotal = numberColumns;
     700  if (inSolve)
     701    numberTotal += numberRows;
     702  currentObj=0.0;
     703  thetaObj=0.0;
     704  for (int iColumn=0;iColumn<numberTotal;iColumn++) {
     705    delta += cost[iColumn]*change[iColumn];
     706    linearCost += cost[iColumn]*solution[iColumn];
     707  }
     708  if (!activated_||!quadraticObjective_) {
     709    currentObj=linearCost;
     710    thetaObj =currentObj + delta*maximumTheta;
     711    if (delta<0.0) {
     712      return maximumTheta;
     713    } else {
     714      printf("odd linear direction %g\n",delta);
     715      return 0.0;
     716    }
     717  }
     718  assert (model);
     719  bool scaling=false;
     720  if ((model->rowScale()||
     721       model->objectiveScale()!=1.0||model->optimizationDirection()!=1.0)&&inSolve)
     722    scaling=true;
     723  const int * columnQuadratic = quadraticObjective_->getIndices();
     724  const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
     725  const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
     726  const double * quadraticElement = quadraticObjective_->getElements();
     727  double a=0.0;
     728  double b=delta;
     729  double c=0.0;
     730  if (!scaling) {
     731    int iColumn;
     732    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     733      double valueI = solution[iColumn];
     734      double changeI = change[iColumn];
     735      CoinBigIndex j;
     736      for (j=columnQuadraticStart[iColumn];
     737           j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     738        int jColumn = columnQuadratic[j];
     739        double valueJ = solution[jColumn];
     740        double changeJ = change[jColumn];
     741        double elementValue = quadraticElement[j];
     742        if (iColumn!=jColumn) {
     743          a += changeI*changeJ*elementValue;
     744          b += (changeI*valueJ+changeJ*valueI)*elementValue;
     745          c += valueI*valueJ*elementValue;
     746        } else {
     747          a += 0.5*changeI*changeI*elementValue;
     748          b += changeI*valueI*elementValue;
     749          c += 0.5*valueI*valueI*elementValue;
     750        }
     751      }
     752    }
     753  } else {
     754    // scaling
     755    const double * columnScale = model->columnScale();
     756    double direction = model->optimizationDirection()*model->objectiveScale();
     757    // direction is actually scale out not scale in
     758    if (direction)
     759      direction = 1.0/direction;
     760    if (!columnScale) {
     761      for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
     762        double valueI = solution[iColumn];
     763        double changeI = change[iColumn];
     764        CoinBigIndex j;
     765        for (j=columnQuadraticStart[iColumn];
     766             j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     767          int jColumn = columnQuadratic[j];
     768          double valueJ = solution[jColumn];
     769          double changeJ = change[jColumn];
     770          double elementValue = quadraticElement[j];
     771          elementValue *= direction;
     772          if (iColumn!=jColumn) {
     773            a += changeI*changeJ*elementValue;
     774            b += (changeI*valueJ+changeJ*valueI)*elementValue;
     775            c += valueI*valueJ*elementValue;
     776          } else {
     777            a += 0.5*changeI*changeI*elementValue;
     778            b += changeI*valueI*elementValue;
     779            c += 0.5*valueI*valueI*elementValue;
     780          }
     781        }
     782      }
     783    } else {
     784      // scaling
     785      for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
     786        double valueI = solution[iColumn];
     787        double changeI = change[iColumn];
     788        double scaleI = columnScale[iColumn]*direction;
     789        CoinBigIndex j;
     790        for (j=columnQuadraticStart[iColumn];
     791             j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     792          int jColumn = columnQuadratic[j];
     793          double valueJ = solution[jColumn];
     794          double changeJ = change[jColumn];
     795          double elementValue = quadraticElement[j];
     796          elementValue *= scaleI*columnScale[jColumn];
     797          if (iColumn!=jColumn) {
     798            a += changeI*changeJ*elementValue;
     799            b += (changeI*valueJ+changeJ*valueI)*elementValue;
     800            c += valueI*valueJ*elementValue;
     801          } else {
     802            a += 0.5*changeI*changeI*elementValue;
     803            b += changeI*valueI*elementValue;
     804            c += 0.5*valueI*valueI*elementValue;
     805          }
     806        }
     807      }
     808    }
     809  }
     810  double theta;
     811  //printf("Current cost %g\n",c+linearCost);
     812  currentObj = c+linearCost;
     813  thetaObj = currentObj + a*maximumTheta*maximumTheta+b*maximumTheta;
     814  // minimize a*x*x + b*x + c
     815  if (a<=0.0) {
     816    theta = maximumTheta;
     817  } else {
     818    theta = -0.5*b/a;
     819  }
     820  predictedObj = currentObj + a*theta*theta+b*theta;
     821  if (b>0.0) {
     822    if (model->messageHandler()->logLevel()&32)
     823      printf("a %g b %g c %g => %g\n",a,b,c,theta);
     824    b=0.0;
     825  }
     826  return min(theta,maximumTheta);
     827}
     828// Scale objective
     829void
     830ClpQuadraticObjective::reallyScale(const double * columnScale)
     831{
     832  const int * columnQuadratic = quadraticObjective_->getIndices();
     833  const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
     834  const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
     835  double * quadraticElement = quadraticObjective_->getMutableElements();
     836  for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
     837    double scaleI = columnScale[iColumn];
     838    objective_[iColumn] *= scaleI;
     839    CoinBigIndex j;
     840    for (j=columnQuadraticStart[iColumn];
     841         j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     842      int jColumn = columnQuadratic[j];
     843      quadraticElement[j] *= scaleI*columnScale[jColumn];
     844    }
     845  }
     846}
     847/* Given a zeroed array sets nonlinear columns to 1.
     848   Returns number of nonlinear columns
     849*/
     850int
     851ClpQuadraticObjective::markNonlinear(char * which)
     852{
     853  int iColumn;
     854  const int * columnQuadratic = quadraticObjective_->getIndices();
     855  const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
     856  const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
     857  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     858    CoinBigIndex j;
     859    for (j=columnQuadraticStart[iColumn];
     860         j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     861      int jColumn = columnQuadratic[j];
     862      which[jColumn]=1;
     863      which[iColumn]=1;
     864    }
     865  }
     866  int numberNonLinearColumns = 0;
     867  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     868    if(which[iColumn])
     869      numberNonLinearColumns++;
     870  }
     871  return numberNonLinearColumns;
     872}
  • trunk/ClpSimplex.cpp

    r389 r393  
    101101  specialOptions_(0),
    102102  lastBadIteration_(-999999),
     103  lastFlaggedIteration_(-999999),
    103104  numberFake_(0),
    104105  progressFlag_(0),
     
    204205  specialOptions_(0),
    205206  lastBadIteration_(-999999),
     207  lastFlaggedIteration_(-999999),
    206208  numberFake_(0),
    207209  progressFlag_(0),
     
    564566ClpSimplex::computeDuals(double * givenDjs)
    565567{
    566   //work space
    567   CoinIndexedVector  * workSpace = rowArray_[0];
    568 
    569   CoinIndexedVector arrayVector;
    570   arrayVector.reserve(numberRows_+1);
    571   CoinIndexedVector previousVector;
    572   previousVector.reserve(numberRows_+1);
    573 
    574 
    575   int iRow;
     568  if (objective_->type()==1||!objective_->activated()) {
     569    // Linear
     570    //work space
     571    CoinIndexedVector  * workSpace = rowArray_[0];
     572   
     573    CoinIndexedVector arrayVector;
     574    arrayVector.reserve(numberRows_+1);
     575    CoinIndexedVector previousVector;
     576    previousVector.reserve(numberRows_+1);
     577   
     578   
     579    int iRow;
    576580#ifdef CLP_DEBUG
    577   workSpace->checkClear();
     581    workSpace->checkClear();
    578582#endif
    579   double * array = arrayVector.denseVector();
    580   int * index = arrayVector.getIndices();
    581   int number=0;
    582   if (!givenDjs) {
    583     for (iRow=0;iRow<numberRows_;iRow++) {
    584       int iPivot=pivotVariable_[iRow];
    585       double value = cost_[iPivot];
    586       if (value) {
    587         array[iRow]=value;
    588         index[number++]=iRow;
    589       }
    590     }
    591   } else {
    592     // dual values pass - djs may not be zero
    593     for (iRow=0;iRow<numberRows_;iRow++) {
    594       int iPivot=pivotVariable_[iRow];
    595       // make sure zero if done
    596       if (!pivoted(iPivot))
    597         givenDjs[iPivot]=0.0;
    598       double value =cost_[iPivot]-givenDjs[iPivot];
    599       if (value) {
    600         array[iRow]=value;
    601         index[number++]=iRow;
    602       }
    603     }
    604   }
    605   arrayVector.setNumElements(number);
    606   // Extended duals before "updateTranspose"
    607   matrix_->dualExpanded(this,&arrayVector,givenDjs,0);
    608 
    609   // Btran basic costs and get as accurate as possible
    610   double lastError=COIN_DBL_MAX;
    611   int iRefine;
    612   double * work = workSpace->denseVector();
    613   CoinIndexedVector * thisVector = &arrayVector;
    614   CoinIndexedVector * lastVector = &previousVector;
    615   factorization_->updateColumnTranspose(workSpace,thisVector);
    616 
    617   for (iRefine=0;iRefine<numberRefinements_+1;iRefine++) {
    618     // check basic reduced costs zero
    619     largestDualError_=0.0;
    620     // would be faster to do just for basic but this reduces code
    621     ClpDisjointCopyN(objectiveWork_,numberColumns_,reducedCostWork_);
    622     transposeTimes(-1.0,array,reducedCostWork_);
    623     // update by duals on sets
    624     matrix_->dualExpanded(this,NULL,NULL,1);
     583    double * array = arrayVector.denseVector();
     584    int * index = arrayVector.getIndices();
     585    int number=0;
    625586    if (!givenDjs) {
    626587      for (iRow=0;iRow<numberRows_;iRow++) {
    627588        int iPivot=pivotVariable_[iRow];
    628         double value;
    629         if (iPivot>=numberColumns_) {
    630           // slack
    631           value = rowObjectiveWork_[iPivot-numberColumns_]
    632             + array[iPivot-numberColumns_];
    633         } else {
    634           // column
    635           value = reducedCostWork_[iPivot];
     589        double value = cost_[iPivot];
     590        if (value) {
     591          array[iRow]=value;
     592          index[number++]=iRow;
    636593        }
    637         work[iRow]=value;
    638         if (fabs(value)>largestDualError_) {
    639           largestDualError_=fabs(value);
    640         }
    641594      }
    642595    } else {
     596      // dual values pass - djs may not be zero
    643597      for (iRow=0;iRow<numberRows_;iRow++) {
    644598        int iPivot=pivotVariable_[iRow];
    645         if (iPivot>=numberColumns_) {
    646           // slack
    647           work[iRow] = rowObjectiveWork_[iPivot-numberColumns_]
    648             + array[iPivot-numberColumns_]-givenDjs[iPivot];
    649         } else {
    650           // column
    651           work[iRow] = reducedCostWork_[iPivot]- givenDjs[iPivot];
    652         }
    653         if (fabs(work[iRow])>largestDualError_) {
    654           largestDualError_=fabs(work[iRow]);
    655           //assert (largestDualError_<1.0e-7);
    656           //if (largestDualError_>1.0e-7)
    657           //printf("large dual error %g\n",largestDualError_);
    658         }
    659       }
    660     }
    661     if (largestDualError_>=lastError) {
    662       // restore
    663       CoinIndexedVector * temp = thisVector;
    664       thisVector = lastVector;
    665       lastVector=temp;
    666       break;
    667     }
    668     if (iRefine<numberRefinements_&&largestDualError_>1.0e-10
    669         &&!givenDjs) {
    670       // try and make better
    671       // save this
    672       CoinIndexedVector * temp = thisVector;
    673       thisVector = lastVector;
    674       lastVector=temp;
    675       int * indexOut = thisVector->getIndices();
    676       int number=0;
    677       array = thisVector->denseVector();
    678       thisVector->clear();
    679       double multiplier = 131072.0;
    680       for (iRow=0;iRow<numberRows_;iRow++) {
    681         double value = multiplier*work[iRow];
     599        // make sure zero if done
     600        if (!pivoted(iPivot))
     601          givenDjs[iPivot]=0.0;
     602        double value =cost_[iPivot]-givenDjs[iPivot];
    682603        if (value) {
    683604          array[iRow]=value;
    684           indexOut[number++]=iRow;
     605          index[number++]=iRow;
     606        }
     607      }
     608    }
     609    arrayVector.setNumElements(number);
     610    // Extended duals before "updateTranspose"
     611    matrix_->dualExpanded(this,&arrayVector,givenDjs,0);
     612   
     613    // Btran basic costs and get as accurate as possible
     614    double lastError=COIN_DBL_MAX;
     615    int iRefine;
     616    double * work = workSpace->denseVector();
     617    CoinIndexedVector * thisVector = &arrayVector;
     618    CoinIndexedVector * lastVector = &previousVector;
     619    factorization_->updateColumnTranspose(workSpace,thisVector);
     620   
     621    for (iRefine=0;iRefine<numberRefinements_+1;iRefine++) {
     622      // check basic reduced costs zero
     623      largestDualError_=0.0;
     624      // would be faster to do just for basic but this reduces code
     625      ClpDisjointCopyN(objectiveWork_,numberColumns_,reducedCostWork_);
     626      transposeTimes(-1.0,array,reducedCostWork_);
     627      // update by duals on sets
     628      matrix_->dualExpanded(this,NULL,NULL,1);
     629      if (!givenDjs) {
     630        for (iRow=0;iRow<numberRows_;iRow++) {
     631          int iPivot=pivotVariable_[iRow];
     632          double value;
     633          if (iPivot>=numberColumns_) {
     634            // slack
     635            value = rowObjectiveWork_[iPivot-numberColumns_]
     636              + array[iPivot-numberColumns_];
     637          } else {
     638            // column
     639            value = reducedCostWork_[iPivot];
     640          }
     641          work[iRow]=value;
     642          if (fabs(value)>largestDualError_) {
     643            largestDualError_=fabs(value);
     644          }
     645        }
     646      } else {
     647        for (iRow=0;iRow<numberRows_;iRow++) {
     648          int iPivot=pivotVariable_[iRow];
     649          if (iPivot>=numberColumns_) {
     650            // slack
     651            work[iRow] = rowObjectiveWork_[iPivot-numberColumns_]
     652              + array[iPivot-numberColumns_]-givenDjs[iPivot];
     653          } else {
     654            // column
     655            work[iRow] = reducedCostWork_[iPivot]- givenDjs[iPivot];
     656          }
     657          if (fabs(work[iRow])>largestDualError_) {
     658            largestDualError_=fabs(work[iRow]);
     659            //assert (largestDualError_<1.0e-7);
     660            //if (largestDualError_>1.0e-7)
     661            //printf("large dual error %g\n",largestDualError_);
     662          }
     663        }
     664      }
     665      if (largestDualError_>=lastError) {
     666        // restore
     667        CoinIndexedVector * temp = thisVector;
     668        thisVector = lastVector;
     669        lastVector=temp;
     670        break;
     671      }
     672      if (iRefine<numberRefinements_&&largestDualError_>1.0e-10
     673          &&!givenDjs) {
     674        // try and make better
     675        // save this
     676        CoinIndexedVector * temp = thisVector;
     677        thisVector = lastVector;
     678        lastVector=temp;
     679        int * indexOut = thisVector->getIndices();
     680        int number=0;
     681        array = thisVector->denseVector();
     682        thisVector->clear();
     683        double multiplier = 131072.0;
     684        for (iRow=0;iRow<numberRows_;iRow++) {
     685          double value = multiplier*work[iRow];
     686          if (value) {
     687            array[iRow]=value;
     688            indexOut[number++]=iRow;
     689            work[iRow]=0.0;
     690          }
    685691          work[iRow]=0.0;
    686692        }
    687         work[iRow]=0.0;
    688       }
    689       thisVector->setNumElements(number);
    690       lastError=largestDualError_;
    691       factorization_->updateColumnTranspose(workSpace,thisVector);
    692       multiplier = 1.0/multiplier;
    693       double * previous = lastVector->denseVector();
    694       number=0;
    695       for (iRow=0;iRow<numberRows_;iRow++) {
    696         double value = previous[iRow] + multiplier*array[iRow];
    697         if (value) {
    698           array[iRow]=value;
    699           indexOut[number++]=iRow;
    700         } else {
    701           array[iRow]=0.0;
     693        thisVector->setNumElements(number);
     694        lastError=largestDualError_;
     695        factorization_->updateColumnTranspose(workSpace,thisVector);
     696        multiplier = 1.0/multiplier;
     697        double * previous = lastVector->denseVector();
     698        number=0;
     699        for (iRow=0;iRow<numberRows_;iRow++) {
     700          double value = previous[iRow] + multiplier*array[iRow];
     701          if (value) {
     702            array[iRow]=value;
     703            indexOut[number++]=iRow;
     704          } else {
     705            array[iRow]=0.0;
     706          }
    702707        }
    703       }
    704       thisVector->setNumElements(number);
    705     } else {
    706       break;
    707     }
    708   }
    709   ClpFillN(work,numberRows_,0.0);
    710   // now look at dual solution
    711   array = thisVector->denseVector();
    712   for (iRow=0;iRow<numberRows_;iRow++) {
    713     // slack
    714     double value = array[iRow];
    715     dual_[iRow]=value;
    716     value += rowObjectiveWork_[iRow];
    717     rowReducedCost_[iRow]=value;
    718   }
    719   ClpDisjointCopyN(objectiveWork_,numberColumns_,reducedCostWork_);
    720   transposeTimes(-1.0,dual_,reducedCostWork_);
    721   // Extended duals and check dual infeasibility
    722   if (!matrix_->skipDualCheck()||algorithm_<0||problemStatus_!=-2)
    723     matrix_->dualExpanded(this,NULL,NULL,2);
    724   // If necessary - override results
    725   if (givenDjs) {
    726     // restore accurate duals
    727     memcpy(givenDjs,dj_,(numberRows_+numberColumns_)*sizeof(double));
    728   }
    729 
     708        thisVector->setNumElements(number);
     709      } else {
     710        break;
     711      }
     712    }
     713    ClpFillN(work,numberRows_,0.0);
     714    // now look at dual solution
     715    array = thisVector->denseVector();
     716    for (iRow=0;iRow<numberRows_;iRow++) {
     717      // slack
     718      double value = array[iRow];
     719      dual_[iRow]=value;
     720      value += rowObjectiveWork_[iRow];
     721      rowReducedCost_[iRow]=value;
     722    }
     723    ClpDisjointCopyN(objectiveWork_,numberColumns_,reducedCostWork_);
     724    transposeTimes(-1.0,dual_,reducedCostWork_);
     725    // Extended duals and check dual infeasibility
     726    if (!matrix_->skipDualCheck()||algorithm_<0||problemStatus_!=-2)
     727      matrix_->dualExpanded(this,NULL,NULL,2);
     728    // If necessary - override results
     729    if (givenDjs) {
     730      // restore accurate duals
     731      memcpy(givenDjs,dj_,(numberRows_+numberColumns_)*sizeof(double));
     732    }
     733  } else {
     734    // Nonlinear
     735    objective_->reducedGradient(this,dj_,false);
     736  }
    730737}
    731738/* Given an existing factorization computes and checks
     
    12661273    solution_[sequenceOut_]=valueOut_;
    12671274  } else {
     1275    //if (objective_->type()<2)
    12681276    assert (fabs(theta_)>1.0e-13);
    12691277    // flip from bound to bound
     
    13021310  int cycle=progress_->cycle(in,out,
    13031311                            directionIn_,directionOut_);
    1304   if (cycle>0) {
     1312  if (cycle>0&&objective_->type()<2) {
     1313    //if (cycle>0) {
    13051314    if (handler_->logLevel()>=63)
    13061315      printf("Cycle of %d\n",cycle);
     
    13471356}
    13481357// Copy constructor.
    1349 ClpSimplex::ClpSimplex(const ClpSimplex &rhs) :
    1350   ClpModel(rhs),
     1358ClpSimplex::ClpSimplex(const ClpSimplex &rhs,int scalingMode) :
     1359  ClpModel(rhs,scalingMode),
    13511360  columnPrimalInfeasibility_(0.0),
    13521361  rowPrimalInfeasibility_(0.0),
     
    14171426  specialOptions_(0),
    14181427  lastBadIteration_(-999999),
     1428  lastFlaggedIteration_(-999999),
    14191429  numberFake_(0),
    14201430  progressFlag_(0),
     
    14441454}
    14451455// Copy constructor from model
    1446 ClpSimplex::ClpSimplex(const ClpModel &rhs) :
    1447   ClpModel(rhs),
     1456ClpSimplex::ClpSimplex(const ClpModel &rhs, int scalingMode) :
     1457  ClpModel(rhs,scalingMode),
    14481458  columnPrimalInfeasibility_(0.0),
    14491459  rowPrimalInfeasibility_(0.0),
     
    15141524  specialOptions_(0),
    15151525  lastBadIteration_(-999999),
     1526  lastFlaggedIteration_(-999999),
    15161527  numberFake_(0),
    15171528  progressFlag_(0),
     
    16551666  specialOptions_ = rhs.specialOptions_;
    16561667  lastBadIteration_ = rhs.lastBadIteration_;
     1668  lastFlaggedIteration_ = rhs.lastFlaggedIteration_;
    16571669  numberFake_ = rhs.numberFake_;
    16581670  progressFlag_ = rhs.progressFlag_;
     
    18411853    }
    18421854  }
     1855  objectiveValue_ += objective_->nonlinearOffset();
    18431856  objectiveValue_ /= (objectiveScale_*rhsScale_);
    18441857}
     
    24802493        rowArray_[iRow]=new CoinIndexedVector();
    24812494        int length =numberRows2+factorization_->maximumPivots();
    2482         if (iRow==3)
     2495        if (iRow==3||objective_->type()>1)
    24832496          length += numberColumns_;
    24842497        rowArray_[iRow]->reserve(length);
     
    32923305int ClpSimplex::dual (int ifValuesPass )
    32933306{
     3307  int saveQuadraticActivated = objective_->activated();
     3308  objective_->setActivated(0);
    32943309  assert (ifValuesPass>=0&&ifValuesPass<3);
    32953310  /*  Note use of "down casting".  The only class the user sees is ClpSimplex.
    3296       Classes ClpSimplexDual, ClpSimplexPrimal, (ClpSimplexPrimalQuadratic)
     3311      Classes ClpSimplexDual, ClpSimplexPrimal, (ClpSimplexNonlinear)
    32973312      and ClpSimplexOther all exist and inherit from ClpSimplex but have no
    32983313      additional data and have no destructor or (non-default) constructor.
     
    33713386      problemStatus_=0;
    33723387  }
     3388  objective_->setActivated(saveQuadraticActivated);
    33733389  return returnCode;
    33743390}
     3391#include "ClpQuadraticObjective.hpp"
    33753392// primal
    33763393int ClpSimplex::primal (int ifValuesPass )
    33773394{
     3395  // See if nonlinear
     3396  if (objective_->type()>1&&objective_->activated())
     3397    return reducedGradient();
     3398 
    33783399  assert (ifValuesPass>=0&&ifValuesPass<3);
    33793400  /*  Note use of "down casting".  The only class the user sees is ClpSimplex.
    3380       Classes ClpSimplexDual, ClpSimplexPrimal, (ClpSimplexPrimalQuadratic)
     3401      Classes ClpSimplexDual, ClpSimplexPrimal, (ClpSimplexNonlinear)
    33813402      and ClpSimplexOther all exist and inherit from ClpSimplex but have no
    33823403      additional data and have no destructor or (non-default) constructor.
     
    35123533  return 0;
    35133534}
    3514 #include "ClpSimplexPrimalQuadratic.hpp"
    3515 /* Solves quadratic problem using SLP - may be used as crash
     3535#include "ClpSimplexNonlinear.hpp"
     3536/* Solves nonlinear problem using SLP - may be used as crash
    35163537   for other algorithms when number of iterations small
    35173538*/
    35183539int
    3519 ClpSimplex::quadraticSLP(int numberPasses, double deltaTolerance)
    3520 {
    3521   return ((ClpSimplexPrimalQuadratic *) this)->primalSLP(numberPasses,deltaTolerance);
    3522 }
    3523 // Solves quadratic using Dantzig's algorithm - primal
    3524 int
    3525 ClpSimplex::quadraticPrimal(int phase)
    3526 {
    3527   return ((ClpSimplexPrimalQuadratic *) this)->primalQuadratic(phase);
     3540ClpSimplex::nonlinearSLP(int numberPasses, double deltaTolerance)
     3541{
     3542  return ((ClpSimplexNonlinear *) this)->primalSLP(numberPasses,deltaTolerance);
     3543}
     3544// Solves non-linear using reduced gradient
     3545int ClpSimplex::reducedGradient(int phase)
     3546{
     3547  if (objective_->type()<2||!objective_->activated()) {
     3548    // no quadratic part
     3549    return primal(0);
     3550  }
     3551  // get feasible
     3552  if ((this->status()<0||numberPrimalInfeasibilities())&&phase==0) {
     3553    objective_->setActivated(0);
     3554    double saveDirection = optimizationDirection();
     3555    setOptimizationDirection(0.0);
     3556    primal(1);
     3557    setOptimizationDirection(saveDirection);
     3558    objective_->setActivated(1);
     3559    // still infeasible
     3560    if (numberPrimalInfeasibilities())
     3561      return 0;
     3562  }
     3563  // Now enter method
     3564  int returnCode = ((ClpSimplexNonlinear *) this)->primal();
     3565  return returnCode;
    35283566}
    35293567#include "ClpPredictorCorrector.hpp"
     
    52605298     
    52615299      if (perturbation_<100) {
    5262         if (algorithm_>0) {
     5300        if (algorithm_>0&&(objective_->type()<2||!objective_->activated())) {
    52635301          ((ClpSimplexPrimal *) this)->perturb(0);
    52645302        } else if (algorithm_<0) {
     
    53645402  status_[sequence] |= 64;
    53655403  matrix_->generalExpanded(this,7,sequence);
     5404  lastFlaggedIteration_=numberIterations_;
    53665405}
    53675406/* Factorizes and returns true if optimal.  Used by user */
  • trunk/ClpSimplexDual.cpp

    r372 r393  
    24582458        progress_->clearBadTimes();
    24592459       
     2460        // Go to safe
     2461        factorization_->pivotTolerance(0.99);
    24602462        forceFactorization_=1; // a bit drastic but ..
    24612463        type = 2;
  • trunk/ClpSimplexPrimal.cpp

    r389 r393  
    621621  int tentativeStatus = problemStatus_;
    622622  int numberThrownOut=1; // to loop round on bad factorization in values pass
     623  double lastSumInfeasibility=COIN_DBL_MAX;
     624  if (numberIterations_)
     625    lastSumInfeasibility=nonLinearCost_->sumInfeasibilities();
    623626  while (numberThrownOut) {
    624627    if (problemStatus_>-3||problemStatus_==-4) {
     
    645648            int saveDense = factorization_->denseThreshold();
    646649            factorization_->setDenseThreshold(0);
     650            // Go to safe
     651            factorization_->pivotTolerance(0.99);
    647652            // make sure will do safe factorization
    648653            pivotVariable_[0]=-1;
     
    663668            forceFactorization_=1; // a bit drastic but ..
    664669            type = 2;
    665             assert (internalFactorize(1)==0);
     670            // Go to safe
     671            factorization_->pivotTolerance(0.99);
     672            if (internalFactorize(1)!=0)
     673               largestPrimalError_=1.0e4; // force other type
    666674          }
    667675          changeMade_++; // say change made
     
    679687    matrix_->generalExpanded(this,9,dummy);
    680688    numberThrownOut=gutsOfSolution(NULL,NULL,(firstFree_>=0));
    681     if (numberThrownOut) {
     689    double sumInfeasibility =  nonLinearCost_->sumInfeasibilities();
     690    if (numberThrownOut||
     691        (sumInfeasibility>1.0e7&&sumInfeasibility>100.0*lastSumInfeasibility
     692         &&factorization_->pivotTolerance()<0.11)) {
    682693      problemStatus_=tentativeStatus;
    683694      doFactorization=true;
     
    10441055                            int valuesPass)
    10451056{
    1046   //rowArray->scanAndPack();
    1047   if (valuesPass) {
     1057  double saveDj = dualIn_;
     1058  if (valuesPass&&objective_->type()<2) {
    10481059    dualIn_ = cost_[sequenceIn_];
    10491060
     
    11511162  //#define CLP_DEBUG
    11521163#ifdef CLP_DEBUG
    1153   if (numberIterations_==3839||numberIterations_==3840) {
     1164  if (numberIterations_==-3839||numberIterations_==-3840) {
    11541165    double dj=cost_[sequenceIn_];
    11551166    printf("cost in on %d is %g, dual in %g\n",sequenceIn_,dj,dualIn_);
     
    12291240
    12301241  bool goBackOne = false;
     1242  if (objective_->type()>1)
     1243    dualIn_=saveDj;
    12311244
    12321245  //printf("%d remain out of %d\n",numberRemaining,number);
     
    13491362                                    -(rhs[iIndex]-spare[iIndex]*theta_));
    13501363      }
    1351 #define CLP_DEBUG
     1364//#define CLP_DEBUG
    13521365#ifdef CLP_DEBUG
    13531366      if (largestInfeasibility>primalTolerance_&&(handler_->logLevel()&32)>-1)
     
    13551368               primalTolerance_,largestInfeasibility);
    13561369#endif
    1357 #undef CLP_DEBUG
     1370//#undef CLP_DEBUG
    13581371      primalTolerance_ = max(primalTolerance_,largestInfeasibility);
    13591372    }
     
    20102023  }
    20112024  numberFlagged += matrix_->generalExpanded(this,8,i);
     2025  if (handler_->logLevel()>2&&numberFlagged&&objective_->type()>1)
     2026    printf("%d unflagged\n",numberFlagged);
    20122027  return numberFlagged;
    20132028}
     
    21152130    if (ifValuesPass) {
    21162131      saveDj=dualIn_;
     2132      //assert (fabs(alpha_)>=1.0e-5||(objective_->type()<2||!objective_->activated())||pivotRow_==-2);
    21172133      if (pivotRow_==-1||(pivotRow_>=0&&fabs(alpha_)<1.0e-5)) {
    2118         if(fabs(dualIn_)<1.0e2*dualTolerance_) {
     2134        if(fabs(dualIn_)<1.0e2*dualTolerance_&&objective_->type()<2) {
    21192135          // try other way
    21202136          directionIn_=-directionIn_;
  • trunk/ClpSolve.cpp

    r390 r393  
    101101  int doCrash=0;
    102102  int doSprint=0;
     103  int doSlp=0;
    103104  if (method!=ClpSolve::useDual&&method!=ClpSolve::useBarrier) {
    104105    switch (options.getSpecialOption(1)) {
     
    152153      doCrash=0;
    153154      doSprint=-1;
     155      break;
     156    case 10:
     157      doIdiot=0;
     158      doCrash=0;
     159      doSprint=0;
     160      if (options.getExtraInfo(1)>0)
     161        doSlp = options.getExtraInfo(1);
    154162      break;
    155163    default:
     
    600608    if (doCrash)
    601609      model2->crash(1000,1);
     610    if (doSlp>0&&objective_->type()==2) {
     611      model2->nonlinearSLP(doSlp,1.0e-5);
     612    }
    602613    model2->primal(1);
    603614    time2 = CoinCpuTime();
     
    12111222      model2->createStatus();
    12121223      model2->dual();
     1224    } else if (maxIts&&quadraticObj) {
     1225      // make sure no status left
     1226      model2->createStatus();
     1227      // solve
     1228      model2->setPerturbation(100);
     1229      model2->reducedGradient(1);
    12131230    }
    12141231    model2->setMaximumIterations(saveMaxIts);
  • trunk/Idiot.cpp

    r257 r393  
    298298#ifndef OSI_IDIOT
    299299  if ((strategy_&32)!=0&&!allOnes) {
    300     scaled = model_->clpMatrix()->scale(model_)==0;
     300    if (model_->scalingFlag()>0)
     301      scaled = model_->clpMatrix()->scale(model_)==0;
    301302    if (scaled) {
    302303      const double * rowScale = model_->rowScale();
  • trunk/Test/ClpMain.cpp

    r389 r393  
    5959
    6060  LOGLEVEL=101,MAXFACTOR,PERTVALUE,MAXITERATION,PRESOLVEPASS,IDIOT,SPRINT,
    61   OUTPUTFORMAT,
     61  OUTPUTFORMAT,SLPVALUE,
    6262 
    6363  DIRECTION=201,DUALPIVOT,SCALING,ERRORSALLOWED,KEEPNAMES,SPARSEFACTOR,
     
    6868  MAXIMIZE,MINIMIZE,EXIT,STDIN,UNITTEST,NETLIB_DUAL,NETLIB_PRIMAL,SOLUTION,
    6969  TIGHTEN,FAKEBOUND,HELP,PLUSMINUS,NETWORK,ALLSLACK,REVERSE,BARRIER,NETLIB_BARRIER,
     70  REALLY_SCALE,
    7071
    7172  INVALID=1000
     
    880881    int doIdiot=-1;
    881882    int outputFormat=2;
     883    int slpValue=-1;
    882884    int doCrash=0;
    883885    int doSprint=-1;
     
    13211323       );
    13221324    parameters[numberParameters++]=
     1325      ClpItem("reallyS!cale","Scales model in place",
     1326              REALLY_SCALE,false);
     1327    parameters[numberParameters++]=
    13231328      ClpItem("restore!Model","Restore model from binary file",
    13241329              RESTORE);
     
    13811386       );
    13821387    parameters[numberParameters-1].setDoubleValue(models->maximumSeconds());
     1388    parameters[numberParameters++]=
     1389      ClpItem("slp!Value","Number of slp passes before primal",
     1390              -1,50000,SLPVALUE,false);
    13831391    parameters[numberParameters++]=
    13841392      ClpItem("sol!ution","Prints solution to file",
     
    15991607            else if (parameters[iParam].type()==OUTPUTFORMAT)
    16001608              outputFormat = value;
     1609            else if (parameters[iParam].type()==SLPVALUE)
     1610              slpValue = value;
    16011611            else
    16021612              parameters[iParam].setIntParameter(models+iModel,value);
     
    17591769              } else if (method==ClpSolve::usePrimalorSprint) {
    17601770                // primal
     1771                // if slp turn everything off
     1772                if (slpValue>0) {
     1773                  doCrash=false;
     1774                  doSprint=0;
     1775                  doIdiot=-1;
     1776                  solveOptions.setSpecialOption(1,10,slpValue); // slp
     1777                  method=ClpSolve::usePrimal;
     1778                }
    17611779                if (doCrash) {
    17621780                  solveOptions.setSpecialOption(1,1); // crash
     
    17661784                } else if (doIdiot>0) {
    17671785                  solveOptions.setSpecialOption(1,2,doIdiot); // idiot
    1768                 } else {
     1786                } else if (slpValue<=0) {
    17691787                  if (doIdiot==0) {
    17701788                    if (doSprint==0)
     
    21742192              assert(dynamic_cast<ClpLinearObjective *> (obj));
    21752193              double offset;
    2176               double * objective = obj->gradient(NULL,offset,true);
     2194              double * objective = obj->gradient(NULL,NULL,offset,true);
    21772195              for (iColumn=0;iColumn<numberColumns;iColumn++) {
    21782196                dualColumnSolution[iColumn] = dualColumnSolution[iColumn];
     
    22852303                  std::endl;
    22862304              }
     2305            }
     2306            break;
     2307          case REALLY_SCALE:
     2308            if (goodModels[iModel]) {
     2309              ClpSimplex newModel(models[iModel],
     2310                                  models[iModel].scalingFlag());
     2311              printf("model really really scaled\n");
     2312              models[iModel]=newModel;
    22872313            }
    22882314            break;
  • trunk/Test/unitTest.cpp

    r389 r393  
    10171017#endif
    10181018  // test network
    1019   //#define QUADRATIC
     1019#define QUADRATIC
    10201020#ifndef QUADRATIC
    10211021  if (1) {   
     
    12261226    CoinRelFltEq eq(1.0e-4);
    12271227    //assert(eq(objValue,820.0));
    1228     solution.setLogLevel(63);
    1229     solution.quadraticPrimal(1);
     1228    //solution.setLogLevel(63);
     1229    solution.primal();
    12301230    int numberRows = solution.numberRows();
    12311231
     
    13041304    int numberColumns=model.numberColumns();
    13051305#if 0
    1306     model.quadraticSLP(50,1.0e-4);
     1306    model.nonlinearSLP(50,1.0e-4);
    13071307#else
    13081308    // Get feasible
     
    13141314    delete saveObjective;
    13151315#endif
    1316     model.setLogLevel(63);
     1316    //model.setLogLevel(63);
    13171317    //exit(77);
    13181318    model.setFactorizationFrequency(10);
    1319     model.quadraticPrimal(1);
     1319    model.primal();
    13201320    double objValue = model.getObjValue();
    1321     const double * solution = model.primalColumnSolution();
    1322     int i;
    1323     for (i=0;i<numberColumns;i++)
    1324       if (solution[i])
    1325         printf("%d %g\n",i,solution[i]);
    13261321    CoinRelFltEq eq(1.0e-4);
    13271322    assert(eq(objValue,-400.92));
     
    13481343    delete [] column;
    13491344    delete [] element;
    1350     solution.quadraticPrimal(1);
    1351     solution.quadraticSLP(50,1.0e-4);
     1345    solution.primal(1);
     1346    solution.nonlinearSLP(50,1.0e-4);
    13521347    double objValue = solution.getObjValue();
    13531348    CoinRelFltEq eq(1.0e-4);
    13541349    assert(eq(objValue,0.5));
    1355     solution.quadraticPrimal();
     1350    solution.primal();
    13561351    objValue = solution.getObjValue();
    13571352    assert(eq(objValue,0.5));
  • trunk/include/ClpInterior.hpp

    r389 r393  
    476476  /// deltaZ.
    477477  double * deltaZ_;
    478   /// deltaT.
    479   double * deltaT_;
     478  /// deltaW.
     479  double * deltaW_;
    480480  /// deltaS.
    481481  double * deltaSU_;
     
    493493  /// rhsZ.
    494494  double * rhsZ_;
    495   /// rhsT.
    496   double * rhsT_;
     495  /// rhsW.
     496  double * rhsW_;
    497497  /// rhs C
    498498  double * rhsC_;
    499499  /// zVec
    500500  double * zVec_;
    501   /// tVec
    502   double * tVec_;
     501  /// wVec
     502  double * wVec_;
    503503  /// cholesky.
    504504  ClpCholeskyBase * cholesky_;
  • trunk/include/ClpLinearObjective.hpp

    r225 r393  
    2222      also returns an offset (to be added to current one)
    2323      If refresh is false then uses last solution
     24      Uses model for scaling
     25      includeLinear 0 - no, 1 as is, 2 as feasible
    2426  */
    25   virtual double * gradient(const double * solution, double & offset,bool refresh);
     27  virtual double * gradient(const ClpSimplex * model,
     28                            const double * solution, double & offset,bool refresh,
     29                            int includeLinear=2);
     30  /** Returns reduced gradient.Returns an offset (to be added to current one).
     31  */
     32  virtual double reducedGradient(ClpSimplex * model, double * region,
     33                                 bool useFeasibleCosts);
     34  /** Returns step length which gives minimum of objective for
     35      solution + theta * change vector up to maximum theta.
     36
     37      arrays are numberColumns+numberRows
     38      Also sets current objective, predicted and at maximumTheta
     39  */
     40  virtual double stepLength(ClpSimplex * model,
     41                            const double * solution,
     42                            const double * change,
     43                            double maximumTheta,
     44                            double & currentObj,
     45                            double & predictedObj,
     46                            double & thetaObj);
    2647  /// Resize objective
    2748  virtual void resize(int newNumberColumns) ;
    2849  /// Delete columns in  objective
    2950  virtual void deleteSome(int numberToDelete, const int * which) ;
     51  /// Scale objective
     52  virtual void reallyScale(const double * columnScale) ;
    3053 
    3154  //@}
  • trunk/include/ClpMatrixBase.hpp

    r389 r393  
    122122  virtual int refresh(ClpSimplex * model)
    123123  { return 0;};
     124  // Really scale matrix
     125  virtual void reallyScale(const double * rowScale, const double * columnScale);
    124126  /** Given positive integer weights for each row fills in sum of weights
    125127      for each column (and slack).
  • trunk/include/ClpMessage.hpp

    r348 r393  
    9595  CLP_BARRIER_STEP,
    9696  CLP_RIM_SCALE,
     97  CLP_SLP_ITER,
    9798  CLP_DUMMY_END
    9899};
  • trunk/include/ClpModel.hpp

    r392 r393  
    4343    ClpModel (  );
    4444
    45     /// Copy constructor.
    46     ClpModel(const ClpModel &);
     45  /** Copy constructor. May scale depending on mode
     46      -1 leave mode as is
     47      0 -off, 1 equilibrium, 2 geometric, 3, auto, 4 dynamic(later)
     48  */
     49    ClpModel(const ClpModel & rhs, int scalingMode=-1);
    4750    /// Assignment operator. This copies the data
    4851    ClpModel & operator=(const ClpModel & rhs);
     
    315318   /// Sets or unsets scaling, 0 -off, 1 equilibrium, 2 geometric, 3, auto, 4 dynamic(later)
    316319   void scaling(int mode=1);
     320  /** If we constructed a "really" scaled model then this reverses the operation.
     321      Quantities may not be exactly as they were before due to rounding errors */
     322  void unscale();
    317323   /// Gets scalingFlag
    318324   inline int scalingFlag() const {return scalingFlag_;};
     
    322328    if (objective_) {
    323329      double offset;
    324       return objective_->gradient(NULL,offset,false);
     330      return objective_->gradient(NULL,NULL,offset,false);
    325331    } else {
    326332      return NULL;
     
    331337    offset=0.0;
    332338    if (objective_) {
    333       return objective_->gradient(solution,offset,refresh);
     339      return objective_->gradient(NULL,solution,offset,refresh);
    334340    } else {
    335341      return NULL;
     
    340346    if (objective_) {
    341347      double offset;
    342       return objective_->gradient(NULL,offset,false);
     348      return objective_->gradient(NULL,NULL,offset,false);
    343349    } else {
    344350      return NULL;
     
    554560                     const double* rowlb, const double* rowub,
    555561                      const double * rowObjective=NULL);
     562  /// Does much of scaling
     563  void gutsOfScaling();
    556564  //@}
    557565
  • trunk/include/ClpObjective.hpp

    r225 r393  
    1212
    1313*/
    14 
     14class ClpSimplex;
     15class ClpModel;
    1516class ClpObjective  {
    1617 
     
    2324      also returns an offset (to be added to current one)
    2425      If refresh is false then uses last solution
     26      Uses model for scaling
     27      includeLinear 0 - no, 1 as is, 2 as feasible
    2528  */
    26   virtual double * gradient(const double * solution, double & offset,bool refresh) = 0;
     29  virtual double * gradient(const ClpSimplex * model,
     30                            const double * solution,
     31                            double & offset,bool refresh,
     32                            int includeLinear=2)=0;
     33  /** Returns reduced gradient.Returns an offset (to be added to current one).
     34  */
     35  virtual double reducedGradient(ClpSimplex * model, double * region,
     36                                 bool useFeasibleCosts)=0;
     37  /** Returns step length which gives minimum of objective for
     38      solution + theta * change vector up to maximum theta.
     39
     40      arrays are numberColumns+numberRows
     41      Also sets current objective, predicted  and at maximumTheta
     42  */
     43  virtual double stepLength(ClpSimplex * model,
     44                            const double * solution,
     45                            const double * change,
     46                            double maximumTheta,
     47                            double & currentObj,
     48                            double & predictedObj,
     49                            double & thetaObj)=0;
    2750  /// Resize objective
    2851  virtual void resize(int newNumberColumns) = 0;
    2952  /// Delete columns in  objective
    3053  virtual void deleteSome(int numberToDelete, const int * which) = 0;
     54  /// Scale objective
     55  virtual void reallyScale(const double * columnScale) =0;
     56  /** Given a zeroed array sets nonlinear columns to 1.
     57      Returns number of nonlinear columns
     58   */
     59  virtual int markNonlinear(char * which);
    3160  //@}
    3261 
     
    6291  inline int type()
    6392  { return type_;};
     93  /// Whether activated
     94  inline int activated() const
     95  {return activated_;};
     96  /// Set whether activated
     97  inline void setActivated(int value)
     98  {activated_=value;};
    6499 
     100  /// Objective offset
     101  inline double nonlinearOffset () const
     102  { return offset_;};
    65103  //@}
    66104
     
    70108  ///@name Protected member data
    71109  //@{
     110  /// Value of non-linear part of objective
     111  double offset_;
    72112  /// Type of objective - linear is 1
    73113  int type_;
     114  /// Whether activated
     115  int activated_;
    74116  //@}
    75117};
  • trunk/include/ClpPackedMatrix.hpp

    r372 r393  
    136136  /// makes sure active columns correct
    137137  virtual int refresh(ClpSimplex * model);
     138  // Really scale matrix
     139  virtual void reallyScale(const double * rowScale, const double * columnScale);
    138140   //@}
    139141
  • trunk/include/ClpPredictorCorrector.hpp

    r389 r393  
    6565                     bool allowIncreasingGap);
    6666  ///:  checks for one step size
    67   bool checkGoodMove2(const double move,double & bestNextGap,
     67  bool checkGoodMove2(double move,double & bestNextGap,
    6868                      bool allowIncreasingGap);
    6969  /// updateSolution.  Updates solution at end of iteration
  • trunk/include/ClpPrimalColumnPivot.hpp

    r256 r393  
    7777  /// Returns true if would not find any column
    7878  virtual bool looksOptimal() const
    79   { return false;};
     79  { return looksOptimal_;};
     80  /// Sets optimality flag (for advanced use)
     81  virtual void setLooksOptimal(bool flag)
     82  { looksOptimal_ = flag;};
    8083  //@}
    8184 
     
    128131  /// Type of column pivot algorithm
    129132  int type_;
     133  /// Says if looks optimal (normally computed)
     134  bool looksOptimal_;
    130135  //@}
    131136};
  • trunk/include/ClpQuadraticObjective.hpp

    r389 r393  
    2323      also returns an offset (to be added to current one)
    2424      If refresh is false then uses last solution
     25      Uses model for scaling
     26      includeLinear 0 - no, 1 as is, 2 as feasible
    2527  */
    26   virtual double * gradient(const double * solution, double & offset,bool refresh);
     28  virtual double * gradient(const ClpSimplex * model,
     29                            const double * solution, double & offset,bool refresh,
     30                            int includeLinear=2);
    2731  /// Resize objective
     32  /** Returns reduced gradient.Returns an offset (to be added to current one).
     33  */
     34  virtual double reducedGradient(ClpSimplex * model, double * region,
     35                                 bool useFeasibleCosts);
     36  /** Returns step length which gives minimum of objective for
     37      solution + theta * change vector up to maximum theta.
     38
     39      arrays are numberColumns+numberRows
     40      Also sets current objective, predicted and at maximumTheta
     41  */
     42  virtual double stepLength(ClpSimplex * model,
     43                            const double * solution,
     44                            const double * change,
     45                            double maximumTheta,
     46                            double & currentObj,
     47                            double & predictedObj,
     48                            double & thetaObj);
    2849  virtual void resize(int newNumberColumns) ;
    2950  /// Delete columns in  objective
    3051  virtual void deleteSome(int numberToDelete, const int * which) ;
     52  /// Scale objective
     53  virtual void reallyScale(const double * columnScale) ;
     54  /** Given a zeroed array sets nonlinear columns to 1.
     55      Returns number of nonlinear columns
     56   */
     57  virtual int markNonlinear(char * which);
    3158 
    3259  //@}
  • trunk/include/ClpSimplex.hpp

    r389 r393  
    7474    ClpSimplex (  );
    7575
    76   /// Copy constructor.
    77   ClpSimplex(const ClpSimplex &);
    78   /// Copy constructor from model.
    79   ClpSimplex(const ClpModel &);
     76  /** Copy constructor. May scale depending on mode
     77      -1 leave mode as is
     78      0 -off, 1 equilibrium, 2 geometric, 3, auto, 4 dynamic(later)
     79  */
     80  ClpSimplex(const ClpSimplex & rhs, int scalingMode =-1);
     81  /** Copy constructor from model. May scale depending on mode
     82      -1 leave mode as is
     83      0 -off, 1 equilibrium, 2 geometric, 3, auto, 4 dynamic(later)
     84  */
     85  ClpSimplex(const ClpModel & rhs, int scalingMode=-1);
    8086  /** Subproblem constructor.  A subset of whole model is created from the
    8187      row and column lists given.  The new order is given by list order and
     
    170176      ifValuesPass==2 just does values pass and then stops */
    171177  int primal(int ifValuesPass=0);
    172   /** Solves quadratic problem using SLP - may be used as crash
     178  /** Solves nonlinear problem using SLP - may be used as crash
    173179      for other algorithms when number of iterations small.
    174180      Also exits if all problematical variables are changing
    175181      less than deltaTolerance
    176182  */
    177   int quadraticSLP(int numberPasses,double deltaTolerance);
    178   /// Solves quadratic using Dantzig's algorithm - primal
    179   int quadraticPrimal(int phase=2);
     183  int nonlinearSLP(int numberPasses,double deltaTolerance);
    180184  /** Solves using barrier (assumes you have good cholesky factor code).
    181185      Does crossover to simplex if asked*/
    182186  int barrier(bool crossover=true);
     187  /** Solves non-linear using reduced gradient.  Phase = 0 get feasible,
     188      =1 use solution */
     189  int reducedGradient(int phase=0);
    183190  /** Dual ranging.
    184191      This computes increase/decrease in cost for each given variable and corresponding
     
    605612  /** Return row or column sections - not as much needed as it
    606613      once was.  These just map into single arrays */
    607   inline double * solutionRegion(int section)
     614  inline double * solutionRegion(int section) const
    608615  { if (!section) return rowActivityWork_; else return columnActivityWork_;};
    609   inline double * djRegion(int section)
     616  inline double * djRegion(int section) const
    610617  { if (!section) return rowReducedCost_; else return reducedCostWork_;};
    611   inline double * lowerRegion(int section)
     618  inline double * lowerRegion(int section) const
    612619  { if (!section) return rowLowerWork_; else return columnLowerWork_;};
    613   inline double * upperRegion(int section)
     620  inline double * upperRegion(int section) const
    614621  { if (!section) return rowUpperWork_; else return columnUpperWork_;};
    615   inline double * costRegion(int section)
     622  inline double * costRegion(int section) const
    616623  { if (!section) return rowObjectiveWork_; else return objectiveWork_;};
    617624  /// Return region as single array
    618   inline double * solutionRegion()
     625  inline double * solutionRegion() const
    619626  { return solution_;};
    620   inline double * djRegion()
     627  inline double * djRegion() const
    621628  { return dj_;};
    622   inline double * lowerRegion()
     629  inline double * lowerRegion() const
    623630  { return lower_;};
    624   inline double * upperRegion()
     631  inline double * upperRegion() const
    625632  { return upper_;};
    626   inline double * costRegion()
     633  inline double * costRegion() const
    627634  { return cost_;};
    628635  inline Status getStatus(int sequence) const
     
    972979  /// So we know when to be cautious
    973980  int lastBadIteration_;
     981  /// So we know when to open up again
     982  int lastFlaggedIteration_;
    974983  /// Can be used for count of fake bounds (dual) or fake costs (primal)
    975984  int numberFake_;
  • trunk/include/ClpSolve.hpp

    r389 r393  
    7474                   8 - do allslack or idiot
    7575                   9 - do allslack or sprint
     76                   10 - slp before
    7677      2 - interrupt handling - 0 yes, 1 no (for threadsafe)
    7778      3 - whether to make +- 1matrix - 0 yes, 1 no
Note: See TracChangeset for help on using the changeset viewer.