Changeset 389


Ignore:
Timestamp:
Jun 15, 2004 1:27:27 PM (16 years ago)
Author:
forrest
Message:

Barrier

Location:
trunk
Files:
22 edited

Legend:

Unmodified
Added
Removed
  • trunk/ClpCholeskyWssmp.cpp

    r369 r389  
    109109// At present I can't get wssmp to work as my libraries seem to be out of sync
    110110// so I have linked in ekkwssmp which is an older version
    111 #if 0
     111//#define WSMP
     112#ifdef WSMP
    112113  extern "C" void wssmp(int * n,
    113114                        int * columnStart , int * rowIndex , double * element,
     
    116117                        double * aux , int * naux ,
    117118                        int   * mrp , int * iparm , double * dparm);
     119extern "C" void wsetmaxthrds(int *);
    118120#else
    119121/* minimum needed for user */
     
    321323  int i0=0;
    322324  int i1=1;
     325#ifdef WSMP   
     326  wsetmaxthrds(&i1);
     327#endif
    323328  wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
    324329         NULL,permuteOut_,permuteIn_,0,&numberRows_,&i1,
     
    413418    }
    414419  }
     420  double delta2 = model_->delta(); // add delta*delta to diagonal
     421  delta2 *= delta2;
    415422  for (iRow=0;iRow<numberRows_;iRow++) {
    416423    double * put = sparseFactor_+choleskyStart_[iRow];
     
    422429      CoinBigIndex startRow=rowStart[iRow];
    423430      CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow];
    424       work[iRow] = diagonalSlack[iRow];
     431      work[iRow] = diagonalSlack[iRow]+delta2;
    425432      for (CoinBigIndex k=startRow;k<endRow;k++) {
    426433        int iColumn=column[k];
     
    455462  //check sizes
    456463  double largest2=maximumAbsElement(sparseFactor_,sizeFactor_);
    457   largest2*=1.0e-19;
     464  largest2*=1.0e-20;
    458465  largest = min (largest2,1.0e-11);
    459466  int numberDroppedBefore=0;
     
    481488  //integerParameters_[11]=1;
    482489  integerParameters_[12]=2;
     490  doubleParameters_[10]=max(1.0e-20,largest2);
     491  if (doubleParameters_[10]>1.0e-3)
     492    integerParameters_[9]=1;
     493  else
     494    integerParameters_[9]=0;
    483495  wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
    484496        NULL,permuteOut_,permuteIn_,NULL,&numberRows_,&i1,
  • trunk/ClpCholeskyWssmpKKT.cpp

    r373 r389  
    1111#include "ClpInterior.hpp"
    1212#include "ClpCholeskyWssmpKKT.hpp"
     13#include "ClpQuadraticObjective.hpp"
    1314#include "ClpMessage.hpp"
    1415
     
    2829    denseThreshold_(denseThreshold)
    2930{
    30   type_=12;
     31  type_=21;
    3132  memset(integerParameters_,0,64*sizeof(int));
    3233  memset(doubleParameters_,0,64*sizeof(double));
     
    8889// At present I can't get wssmp to work as my libraries seem to be out of sync
    8990// so I have linked in ekkwssmp which is an older version
    90 #if 0
     91//#define WSMP
     92#ifdef WSMP
    9193  extern "C" void wssmp(int * n,
    9294                        int * columnStart , int * rowIndex , double * element,
     
    9597                        double * aux , int * naux ,
    9698                        int   * mrp , int * iparm , double * dparm);
     99extern "C" void wsetmaxthrds(int *);
    97100#else
    98101/* minimum needed for user */
     
    145148  numberRowsDropped_=0;
    146149  model_=model;
     150  CoinPackedMatrix * quadratic = NULL;
     151  ClpQuadraticObjective * quadraticObj =
     152    (dynamic_cast< ClpQuadraticObjective*>(model_->objectiveAsObject()));
     153  if (quadraticObj)
     154    quadratic = quadraticObj->quadraticObjective();
    147155  int numberElements = model_->clpMatrix()->getNumElements();
    148156  numberElements = numberElements+2*numberRowsModel+numberTotal;
     157  if (quadratic)
     158    numberElements += quadratic->getNumElements();
    149159  // Space for starts
    150160  choleskyStart_ = new CoinBigIndex[numberRows_+1];
     
    178188  sizeFactor_=0;
    179189  // matrix
    180   for (iColumn=0;iColumn<numberColumns;iColumn++) {
    181     choleskyStart_[iColumn]=sizeFactor_;
    182     choleskyRow_[sizeFactor_++]=iColumn;
    183     CoinBigIndex start=columnStart[iColumn];
    184     CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
    185     for (CoinBigIndex j=start;j<end;j++) {
    186       choleskyRow_[sizeFactor_++]=row[j]+numberTotal;
     190  if (!quadratic) {
     191    for (iColumn=0;iColumn<numberColumns;iColumn++) {
     192      choleskyStart_[iColumn]=sizeFactor_;
     193      choleskyRow_[sizeFactor_++]=iColumn;
     194      CoinBigIndex start=columnStart[iColumn];
     195      CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
     196      for (CoinBigIndex j=start;j<end;j++) {
     197        choleskyRow_[sizeFactor_++]=row[j]+numberTotal;
     198      }
     199    }
     200  } else {
     201    // Quadratic
     202    const int * columnQuadratic = quadratic->getIndices();
     203    const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
     204    const int * columnQuadraticLength = quadratic->getVectorLengths();
     205    //const double * quadraticElement = quadratic->getElements();
     206    for (iColumn=0;iColumn<numberColumns;iColumn++) {
     207      choleskyStart_[iColumn]=sizeFactor_;
     208      choleskyRow_[sizeFactor_++]=iColumn;
     209      for (CoinBigIndex j=columnQuadraticStart[iColumn];
     210           j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     211        int jColumn = columnQuadratic[j];
     212        if (jColumn>iColumn)
     213          choleskyRow_[sizeFactor_++]=jColumn;
     214      }
     215      CoinBigIndex start=columnStart[iColumn];
     216      CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
     217      for (CoinBigIndex j=start;j<end;j++) {
     218        choleskyRow_[sizeFactor_++]=row[j]+numberTotal;
     219      }
    187220    }
    188221  }
     
    193226    choleskyRow_[sizeFactor_++]=iColumn-numberColumns+numberTotal;
    194227  }
    195   // Transpose - zero diagonal (may regularize)
     228  // Transpose - nonzero diagonal (may regularize)
    196229  for (iRow=0;iRow<numberRowsModel;iRow++) {
    197230    choleskyStart_[iRow+numberTotal]=sizeFactor_;
     
    205238  int i0=0;
    206239  int i1=1;
     240#ifdef WSMP   
     241  wsetmaxthrds(&i1);
     242#endif
    207243  wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
    208244         NULL,permuteOut_,permuteIn_,0,&numberRows_,&i1,
     
    279315 
    280316  CoinBigIndex numberElements=0;
     317  CoinPackedMatrix * quadratic = NULL;
     318  ClpQuadraticObjective * quadraticObj =
     319    (dynamic_cast< ClpQuadraticObjective*>(model_->objectiveAsObject()));
     320  if (quadraticObj)
     321    quadratic = quadraticObj->quadraticObjective();
    281322  // matrix
     323  largest=1.0e-100;
     324  if (!quadratic) {
     325    for (iColumn=0;iColumn<numberColumns;iColumn++) {
     326      choleskyStart_[iColumn]=numberElements;
     327      double value = diagonal[iColumn];
     328      if (fabs(value)>1.0e-100) {
     329        value = 1.0/value;
     330        largest = max(largest,fabs(value));
     331        sparseFactor_[numberElements] = -value;
     332        choleskyRow_[numberElements++]=iColumn;
     333        CoinBigIndex start=columnStart[iColumn];
     334        CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
     335        for (CoinBigIndex j=start;j<end;j++) {
     336          choleskyRow_[numberElements]=row[j]+numberTotal;
     337          sparseFactor_[numberElements++]=element[j];
     338        }
     339      } else {
     340        sparseFactor_[numberElements] = -1.0e100;
     341        choleskyRow_[numberElements++]=iColumn;
     342      }
     343    }
     344  } else {
     345    // Quadratic
     346    const int * columnQuadratic = quadratic->getIndices();
     347    const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
     348    const int * columnQuadraticLength = quadratic->getVectorLengths();
     349    const double * quadraticElement = quadratic->getElements();
     350    for (iColumn=0;iColumn<numberColumns;iColumn++) {
     351      choleskyStart_[iColumn]=numberElements;
     352      CoinBigIndex savePosition = numberElements;
     353      choleskyRow_[numberElements++]=iColumn;
     354      double value = diagonal[iColumn];
     355      if (fabs(value)>1.0e-100) {
     356        value = 1.0/value;
     357        for (CoinBigIndex j=columnQuadraticStart[iColumn];
     358             j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     359          int jColumn = columnQuadratic[j];
     360          if (jColumn>iColumn) {
     361            sparseFactor_[numberElements]=quadraticElement[j];
     362            choleskyRow_[numberElements++]=jColumn;
     363          } else if (iColumn==jColumn) {
     364            value += quadraticElement[j];
     365          }
     366        }
     367        largest = max(largest,fabs(value));
     368        sparseFactor_[savePosition] = -value;
     369        CoinBigIndex start=columnStart[iColumn];
     370        CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
     371        for (CoinBigIndex j=start;j<end;j++) {
     372          choleskyRow_[numberElements]=row[j]+numberTotal;
     373          sparseFactor_[numberElements++]=element[j];
     374        }
     375      } else {
     376        value = 1.0e100;
     377        sparseFactor_[savePosition] = -value;
     378      }
     379    }
     380  }
    282381  for (iColumn=0;iColumn<numberColumns;iColumn++) {
     382    assert (sparseFactor_[choleskyStart_[iColumn]]<0.0);
     383  }
     384  // slacks
     385  for (iColumn=numberColumns;iColumn<numberTotal;iColumn++) {
    283386    choleskyStart_[iColumn]=numberElements;
    284387    double value = diagonal[iColumn];
    285     if (fabs(value)>1.0e-100)
     388    if (fabs(value)>1.0e-100) {
    286389      value = 1.0/value;
    287     else
     390      largest = max(largest,fabs(value));
     391    } else {
    288392      value = 1.0e100;
    289     sparseFactor_[numberElements] = -value;
    290     choleskyRow_[numberElements++]=iColumn;
    291     CoinBigIndex start=columnStart[iColumn];
    292     CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
    293     for (CoinBigIndex j=start;j<end;j++) {
    294       choleskyRow_[numberElements]=row[j]+numberTotal;
    295       sparseFactor_[numberElements++]=element[j];
    296     }
    297   }
    298   // slacks
    299   for (;iColumn<numberTotal;iColumn++) {
    300     choleskyStart_[iColumn]=numberElements;
    301     double value = diagonal[iColumn];
    302     if (fabs(value)>1.0e-100)
    303       value = 1.0/value;
    304     else
    305       value = 1.0e100;
     393    }
    306394    sparseFactor_[numberElements] = -value;
    307395    choleskyRow_[numberElements++]=iColumn;
     
    310398  }
    311399  // Finish diagonal
     400  double delta2 = model_->delta(); // add delta*delta to bottom
     401  delta2 *= delta2;
    312402  for (iRow=0;iRow<numberRowsModel;iRow++) {
    313403    choleskyStart_[iRow+numberTotal]=numberElements;
    314404    choleskyRow_[numberElements]=iRow+numberTotal;
    315     sparseFactor_[numberElements++]=0.0;
     405    sparseFactor_[numberElements++]=delta2;
    316406  }
    317407  choleskyStart_[numberRows_]=numberElements;
     
    326416  integerParameters_[30]=1;
    327417  doubleParameters_[20]=1.0e100;
     418#ifndef WSMP
     419  // Set up LDL cutoff
     420  integerParameters_[34]=numberTotal;
     421  doubleParameters_[10]=min(largest*1.0e-19,1.0e-14);
     422  doubleParameters_[20]=1.0e-15;
     423  doubleParameters_[34]=1.0e-12;
     424  //printf("tol is %g\n",doubleParameters_[10]);
     425  //doubleParameters_[10]=1.0e-17;
     426#endif
    328427  int * rowsDropped2 = new int[numberRows_];
    329428  CoinZeroN(rowsDropped2,numberRows_);
     
    337436  newDropped=integerParameters_[20];
    338437#if 1
     438  // Should save adjustments in ..R_
    339439  int n1=0,n2=0;
     440  double * primalR = model_->primalR();
     441  double * dualR = model_->dualR();
    340442  for (iRow=0;iRow<numberTotal;iRow++) {
    341443    if (rowsDropped2[iRow]) {
    342444      n1++;
    343445      //printf("row region1 %d dropped\n",iRow);
    344       rowsDropped_[iRow]=1;
     446      //rowsDropped_[iRow]=1;
     447      rowsDropped_[iRow]=0;
     448      primalR[iRow]=doubleParameters_[20];
    345449    } else {
    346450      rowsDropped_[iRow]=0;
     451      primalR[iRow]=0.0;
    347452    }
    348453  }
     
    351456      n2++;
    352457      //printf("row region2 %d dropped\n",iRow);
    353       rowsDropped_[iRow]=1;
     458      //rowsDropped_[iRow]=1;
     459      rowsDropped_[iRow]=0;
     460      dualR[iRow-numberTotal]=doubleParameters_[34];
    354461    } else {
    355462      rowsDropped_[iRow]=0;
    356     }
    357   }
    358   printf("%d rows dropped in region1, %d in region2\n",n1,n2);
     463      dualR[iRow-numberTotal]=0.0;
     464    }
     465  }
     466  //printf("%d rows dropped in region1, %d in region2\n",n1,n2);
    359467#endif
    360468  delete [] rowsDropped2;
     
    414522#endif
    415523  CoinMemcpyN(array+numberTotal,numberRowsModel,region2);
    416 #if 0
     524#if 1
    417525  CoinMemcpyN(array,numberTotal,region1);
    418526#else
  • trunk/ClpInterior.cpp

    r388 r389  
    2121#include "ClpCholeskyDense.hpp"
    2222#include "ClpLinearObjective.hpp"
     23#include "ClpQuadraticObjective.hpp"
    2324#include <cfloat>
    2425
     
    9091  deltaSU_(NULL),
    9192  deltaSL_(NULL),
     93  primalR_(NULL),
     94  dualR_(NULL),
    9295  rhsB_(NULL),
    9396  rhsU_(NULL),
     
    177180    deltaSU_(NULL),
    178181    deltaSL_(NULL),
     182    primalR_(NULL),
     183    dualR_(NULL),
    179184    rhsB_(NULL),
    180185    rhsU_(NULL),
     
    250255  deltaSU_(NULL),
    251256  deltaSL_(NULL),
     257  primalR_(NULL),
     258  dualR_(NULL),
    252259  rhsB_(NULL),
    253260  rhsU_(NULL),
     
    326333  deltaSU_(NULL),
    327334  deltaSL_(NULL),
     335  primalR_(NULL),
     336  dualR_(NULL),
    328337  rhsB_(NULL),
    329338  rhsU_(NULL),
     
    421430  deltaSU_ = ClpCopyOfArray(rhs.deltaSU_,numberRows_+numberColumns_);
    422431  deltaSL_ = ClpCopyOfArray(rhs.deltaSL_,numberRows_+numberColumns_);
     432  primalR_ = ClpCopyOfArray(rhs.primalR_,numberRows_+numberColumns_);
     433  dualR_ = ClpCopyOfArray(rhs.dualR_,numberRows_+numberColumns_);
    423434  rhsB_ = ClpCopyOfArray(rhs.rhsB_,numberRows_);
    424435  rhsU_ = ClpCopyOfArray(rhs.rhsU_,numberRows_+numberColumns_);
     
    489500  delete [] deltaSL_;
    490501  deltaSL_ = NULL;
     502  delete [] primalR_;
     503  primalR_=NULL;
     504  delete [] dualR_;
     505  dualR_=NULL;
    491506  delete [] rhsB_;
    492507  rhsB_ = NULL;
     
    530545  cost_ = new double[nTotal];
    531546  int i;
    532   double direction = optimizationDirection_;
     547  double direction = optimizationDirection_*objectiveScale_;
    533548  // direction is actually scale out not scale in
    534549  if (direction)
     
    539554  memset(cost_+numberColumns_,0,numberRows_*sizeof(double));
    540555  // do scaling if needed
    541   if (scalingFlag_>0&&!rowScale_&&0) {
     556  if (scalingFlag_>0&&!rowScale_) {
    542557    if (matrix_->scale(this))
    543558      scalingFlag_=-scalingFlag_; // not scaled after all
     
    572587  if (!sanityCheck())
    573588    goodMatrix=false;
    574   if (rowScale_) {
    575     for (i=0;i<numberColumns_;i++)
     589  if(rowScale_) {
     590    for (i=0;i<numberColumns_;i++) {
     591      double multiplier = rhsScale_/columnScale_[i];
    576592      cost_[i] *= columnScale_[i];
    577     for (i=0;i<numberColumns_;i++) {
    578       double multiplier = 1.0/columnScale_[i];
    579593      if (columnLowerWork_[i]>-1.0e50)
    580594        columnLowerWork_[i] *= multiplier;
     
    584598    }
    585599    for (i=0;i<numberRows_;i++) {
    586       double multiplier = rowScale_[i];
     600      double multiplier = rhsScale_*rowScale_[i];
    587601      if (rowLowerWork_[i]>-1.0e50)
    588602        rowLowerWork_[i] *= multiplier;
    589603      if (rowUpperWork_[i]<1.0e50)
    590604        rowUpperWork_[i] *= multiplier;
     605    }
     606  } else if (rhsScale_!=1.0) {
     607    for (i=0;i<numberColumns_+numberRows_;i++) {
     608      if (lower_[i]>-1.0e50)
     609        lower_[i] *= rhsScale_;
     610      if (upper_[i]<1.0e50)
     611        upper_[i] *= rhsScale_;
     612     
    591613    }
    592614  }
     
    619641  deltaSL_ = new double [nTotal];
    620642  CoinZeroN(deltaSL_,nTotal);
     643  assert (!primalR_);
     644  assert (!dualR_);
     645  // create arrays if we are doing KKT
     646  if (cholesky_->type()>20) {
     647    primalR_ = new double [nTotal];
     648    CoinZeroN(primalR_,nTotal);
     649    dualR_ = new double [numberRows_];
     650    CoinZeroN(dualR_,numberRows_);
     651  }
    621652  assert (!rhsB_);
    622653  rhsB_ = new double [numberRows_];
     
    657688{
    658689  int i;
    659   if (optimizationDirection_!=1.0) {
     690  if (optimizationDirection_!=1.0||objectiveScale_!=1.0) {
     691    double scaleC = optimizationDirection_/objectiveScale_;
    660692    // and modify all dual signs
    661693    for (i=0;i<numberColumns_;i++)
    662       reducedCost_[i] = optimizationDirection_*dj_[i];
     694      reducedCost_[i] = scaleC*dj_[i];
    663695    for (i=0;i<numberRows_;i++)
    664       dual_[i] *= optimizationDirection_;
     696      dual_[i] *= scaleC;
    665697  }
    666698  if (rowScale_) {
     699    double scaleR = 1.0/rhsScale_;
    667700    for (i=0;i<numberColumns_;i++) {
    668701      double scaleFactor = columnScale_[i];
    669       columnActivity_[i] *= scaleFactor;
    670       reducedCost_[i] /= scaleFactor;
     702      double valueScaled = columnActivity_[i];
     703      columnActivity_[i] = valueScaled*scaleFactor*scaleR;
     704      double valueScaledDual = reducedCost_[i];
     705      reducedCost_[i] = valueScaledDual/scaleFactor;
    671706    }
    672707    for (i=0;i<numberRows_;i++) {
    673708      double scaleFactor = rowScale_[i];
    674       rowActivity_[i] /= scaleFactor;
    675       dual_[i] *= scaleFactor;
     709      double valueScaled = rowActivity_[i];
     710      rowActivity_[i] = (valueScaled*scaleR)/scaleFactor;
     711      double valueScaledDual = dual_[i];
     712      dual_[i] = valueScaledDual*scaleFactor;
     713    }
     714  } else if (rhsScale_!=1.0) {
     715    double scaleR = 1.0/rhsScale_;
     716    for (i=0;i<numberColumns_;i++) {
     717      double valueScaled = columnActivity_[i];
     718      columnActivity_[i] = valueScaled*scaleR;
     719    }
     720    for (i=0;i<numberRows_;i++) {
     721      double valueScaled = rowActivity_[i];
     722      rowActivity_[i] = valueScaled*scaleR;
    676723    }
    677724  }
     
    946993{
    947994  int iRow,iColumn;
     995  CoinMemcpyN(cost_,numberColumns_,reducedCost_);
     996  matrix_->transposeTimes(-1.0,dual_,reducedCost_);
     997  // Now modify reduced costs for quadratic
     998  double quadraticOffset=quadraticDjs(reducedCost_,
     999                                      solution_,scaleFactor_);
    9481000
    9491001  objectiveValue_ = 0.0;
     
    9531005  double dualTolerance =  10.0*dblParam_[ClpDualTolerance];
    9541006  double primalTolerance =  dblParam_[ClpPrimalTolerance];
     1007  double primalTolerance2 =  10.0*dblParam_[ClpPrimalTolerance];
    9551008  worstComplementarity_=0.0;
    9561009  complementarityGap_=0.0;
    9571010
     1011  // Done scaled - use permanent regions for output
     1012  // but internal for bounds
     1013  const double * lower = lower_+numberColumns_;
     1014  const double * upper = upper_+numberColumns_;
    9581015  for (iRow=0;iRow<numberRows_;iRow++) {
    9591016    double infeasibility=0.0;
    960     double distanceUp = min(rowUpper_[iRow]-
     1017    double distanceUp = min(upper[iRow]-
    9611018      rowActivity_[iRow],1.0e10);
    9621019    double distanceDown = min(rowActivity_[iRow] -
    963       rowLower_[iRow],1.0e10);
    964     if (distanceUp>primalTolerance) {
     1020      lower[iRow],1.0e10);
     1021    if (distanceUp>primalTolerance2) {
    9651022      double value = dual_[iRow];
    9661023      // should not be negative
    9671024      if (value<-dualTolerance) {
     1025        sumDualInfeasibilities_ += -dualTolerance-value;
    9681026        value = - value*distanceUp;
    9691027        if (value>worstComplementarity_)
     
    9721030      }
    9731031    }
    974     if (distanceDown>primalTolerance) {
     1032    if (distanceDown>primalTolerance2) {
    9751033      double value = dual_[iRow];
    9761034      // should not be positive
    9771035      if (value>dualTolerance) {
     1036        sumDualInfeasibilities_ += value-dualTolerance;
    9781037        value =  value*distanceDown;
    9791038        if (value>worstComplementarity_)
     
    9821041      }
    9831042    }
    984     if (rowActivity_[iRow]>rowUpper_[iRow]) {
    985       infeasibility=rowActivity_[iRow]-rowUpper_[iRow];
    986     } else if (rowActivity_[iRow]<rowLower_[iRow]) {
    987       infeasibility=rowLower_[iRow]-rowActivity_[iRow];
     1043    if (rowActivity_[iRow]>upper[iRow]) {
     1044      infeasibility=rowActivity_[iRow]-upper[iRow];
     1045    } else if (rowActivity_[iRow]<lower[iRow]) {
     1046      infeasibility=lower[iRow]-rowActivity_[iRow];
    9881047    }
    9891048    if (infeasibility>primalTolerance) {
     
    9911050    }
    9921051  }
     1052  lower = lower_;
     1053  upper = upper_;
    9931054  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    9941055    double infeasibility=0.0;
    9951056    objectiveValue_ += cost_[iColumn]*columnActivity_[iColumn];
    996     double distanceUp = min(columnUpper_[iColumn]-
     1057    double distanceUp = min(upper[iColumn]-
    9971058      columnActivity_[iColumn],1.0e10);
    9981059    double distanceDown = min(columnActivity_[iColumn] -
    999       columnLower_[iColumn],1.0e10);
    1000     if (distanceUp>primalTolerance) {
     1060      lower[iColumn],1.0e10);
     1061    if (distanceUp>primalTolerance2) {
    10011062      double value = reducedCost_[iColumn];
    10021063      // should not be negative
    10031064      if (value<-dualTolerance) {
     1065        sumDualInfeasibilities_ += -dualTolerance-value;
    10041066        value = - value*distanceUp;
    10051067        if (value>worstComplementarity_)
     
    10081070      }
    10091071    }
    1010     if (distanceDown>primalTolerance) {
     1072    if (distanceDown>primalTolerance2) {
    10111073      double value = reducedCost_[iColumn];
    10121074      // should not be positive
    10131075      if (value>dualTolerance) {
     1076        sumDualInfeasibilities_ += value-dualTolerance;
    10141077        value =  value*distanceDown;
    10151078        if (value>worstComplementarity_)
     
    10181081      }
    10191082    }
    1020     if (columnActivity_[iColumn]>columnUpper_[iColumn]) {
    1021       infeasibility=columnActivity_[iColumn]-columnUpper_[iColumn];
    1022     } else if (columnActivity_[iColumn]<columnLower_[iColumn]) {
    1023       infeasibility=columnLower_[iColumn]-columnActivity_[iColumn];
     1083    if (columnActivity_[iColumn]>upper[iColumn]) {
     1084      infeasibility=columnActivity_[iColumn]-upper[iColumn];
     1085    } else if (columnActivity_[iColumn]<lower[iColumn]) {
     1086      infeasibility=lower[iColumn]-columnActivity_[iColumn];
    10241087    }
    10251088    if (infeasibility>primalTolerance) {
     
    10271090    }
    10281091  }
     1092  // add in offset
     1093  objectiveValue_ += 0.5*quadraticOffset;
    10291094}
    10301095// Set cholesky (and delete present one)
     
    10761141// fix variables interior says should be
    10771142void
    1078 ClpInterior::fixFixed()
     1143ClpInterior::fixFixed(bool reallyFix)
    10791144{
    10801145  int i;
     
    10841149        if (fixedOrFree(i)) {
    10851150          if (columnActivity_[i]-columnLower_[i]<columnUpper_[i]-columnActivity_[i]) {
    1086             columnUpper_[i]=columnLower_[i];
     1151            if (reallyFix)
     1152              columnUpper_[i]=columnLower_[i];
    10871153            columnActivity_[i]=columnLower_[i];
    10881154          } else {
    1089             columnLower_[i]=columnUpper_[i];
     1155            if (reallyFix)
     1156              columnLower_[i]=columnUpper_[i];
    10901157            columnActivity_[i]=columnUpper_[i];
    10911158          }
     
    11101177  }
    11111178}
     1179/* Modifies djs to allow for quadratic.
     1180   returns quadratic offset */
     1181double
     1182ClpInterior::quadraticDjs(double * djRegion, const double * solution,
     1183                          double scaleFactor)
     1184{
     1185  double quadraticOffset=0.0;
     1186#ifndef NO_RTTI
     1187  ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
     1188#else
     1189  ClpQuadraticObjective * quadraticObj = NULL;
     1190  if (objective_->type()==2)
     1191    quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
     1192#endif
     1193  if (quadraticObj) {
     1194    CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
     1195    const int * columnQuadratic = quadratic->getIndices();
     1196    const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
     1197    const int * columnQuadraticLength = quadratic->getVectorLengths();
     1198    double * quadraticElement = quadratic->getMutableElements();
     1199    int numberColumns = quadratic->getNumCols();
     1200    for (int iColumn=0;iColumn<numberColumns;iColumn++) {
     1201      double value = 0.0;
     1202      for (CoinBigIndex j=columnQuadraticStart[iColumn];
     1203           j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     1204        int jColumn = columnQuadratic[j];
     1205        double valueJ = solution[jColumn];
     1206        double elementValue = quadraticElement[j];
     1207        //value += valueI*valueJ*elementValue;
     1208        value += valueJ*elementValue;
     1209        quadraticOffset += solution[iColumn]*valueJ*elementValue;
     1210      }
     1211      djRegion[iColumn] += scaleFactor*value;
     1212    }
     1213  }
     1214  return quadraticOffset;
     1215}
  • trunk/ClpModel.cpp

    r381 r389  
    3131  objectiveValue_(0.0),
    3232  smallElement_(1.0e-20),
     33  objectiveScale_(1.0),
     34  rhsScale_(1.0),
    3335  numberRows_(0),
    3436  numberColumns_(0),
     
    336338  objectiveValue_=rhs.objectiveValue_;
    337339  smallElement_ = rhs.smallElement_;
     340  objectiveScale_ = rhs.objectiveScale_;
     341  rhsScale_ = rhs.rhsScale_;
    338342  numberIterations_ = rhs.numberIterations_;
    339343  solveType_ = rhs.solveType_;
     
    11991203      integerType_ = NULL;
    12001204    }
     1205    // get quadratic part
     1206    if (m.reader()->whichSection (  ) == COIN_QUAD_SECTION ) {
     1207      int * start=NULL;
     1208      int * column = NULL;
     1209      double * element = NULL;
     1210      status=m.readQuadraticMps(NULL,start,column,element,2);
     1211      if (!status||ignoreErrors)
     1212        loadQuadraticObjective(numberColumns_,start,column,element);
     1213      delete [] start;
     1214      delete [] column;
     1215      delete [] element;
     1216    }
    12011217    // set problem name
    12021218    setStrParam(ClpProbName,m.getProblemName());
     
    14361452  objectiveValue_=rhs->objectiveValue_;
    14371453  smallElement_ = rhs->smallElement_;
     1454  objectiveScale_ = rhs->objectiveScale_;
     1455  rhsScale_ = rhs->rhsScale_;
    14381456  numberIterations_ = rhs->numberIterations_;
    14391457  solveType_ = rhs->solveType_;
  • trunk/ClpPredictorCorrector.cpp

    r382 r389  
    1717#include "ClpCholeskyBase.hpp"
    1818#include "ClpHelperFunctions.hpp"
     19#include "ClpQuadraticObjective.hpp"
    1920#include <cfloat>
    2021#include <cassert>
     
    4445  }
    4546  ClpMatrixBase * saveMatrix = NULL;
     47  // If quadratic then make copy so we can actually scale or normalize
     48#ifndef NO_RTTI
     49  ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
     50#else
     51  ClpQuadraticObjective * quadraticObj = NULL;
     52  if (objective_->type()==2)
     53    quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
     54#endif
     55  /* If modeSwitch is :
     56     0 - normal
     57     1 - bit switch off centering
     58     2 - bit always do type 2
     59     4 - accept corrector nearly always
     60  */
     61  int modeSwitch=0;
     62  //if (quadraticObj)
     63  //modeSwitch |= 1; // switch off centring for now
     64  //if (quadraticObj)
     65  //modeSwitch |=4;
     66  ClpObjective * saveObjective = NULL;
     67  if (quadraticObj) {
     68    saveObjective=objective_;
     69    // We are going to make matrix full rather half
     70    objective_ = new ClpQuadraticObjective(*quadraticObj,1);
     71  }
     72  bool allowIncreasingGap = (modeSwitch&4)!=0;
    4673  // If scaled then really scale matrix
    4774  if (scalingFlag_>0&&rowScale_) {
     
    112139  const double tau   = 0.00002;
    113140  double lastComplementarityGap=COIN_DBL_MAX;
     141  double lastStep=1.0;
    114142  // use to see if to take affine
    115143  double checkGap = COIN_DBL_MAX;
     
    125153  double * saveZ = new double[numberTotal];
    126154  double * saveT = new double[numberTotal];
     155  double * saveSL = new double[numberTotal];
     156  double * saveSU = new double[numberTotal];
    127157  // Save smallest mu used in primal dual moves
    128158  double smallestPrimalDualMu=COIN_DBL_MAX;
     159  double objScale = optimizationDirection_/
     160    (rhsScale_*objectiveScale_);
    129161  while (problemStatus_<0) {
    130162    //#define FULL_DEBUG
     
    145177    complementarityGap_=complementarityGap(numberComplementarityPairs_,
    146178                                           numberComplementarityItems_,0);
    147       handler_->message(CLP_BARRIER_ITERATION,messages_)
    148     <<numberIterations_
    149     <<primalObjective_*optimizationDirection_- dblParam_[ClpObjOffset]
    150     << dualObjective_*optimizationDirection_- dblParam_[ClpObjOffset]
    151     <<complementarityGap_
    152     <<numberFixedTotal
    153     <<cholesky_->rank()
    154     <<CoinMessageEol;
     179    handler_->message(CLP_BARRIER_ITERATION,messages_)
     180      <<numberIterations_
     181      <<primalObjective_*objScale- dblParam_[ClpObjOffset]
     182      << dualObjective_*objScale- dblParam_[ClpObjOffset]
     183      <<complementarityGap_
     184      <<numberFixedTotal
     185      <<cholesky_->rank()
     186      <<CoinMessageEol;
     187    // move up history
     188    for (i=1;i<LENGTH_HISTORY;i++)
     189      historyInfeasibility_[i-1]=historyInfeasibility_[i];
     190    historyInfeasibility_[LENGTH_HISTORY-1]=complementarityGap_;
     191    // switch off saved if changes
     192    if (saveIteration+10<numberIterations_&&
     193        complementarityGap_*2.0<historyInfeasibility_[0])
     194      saveIteration=-1;
     195    lastStep = min(actualPrimalStep_,actualDualStep_);
    155196    double goodGapChange;
    156197    if (!sloppyOptimal) {
     
    216257          if (maximumBoundInfeasibility_>1.0e-1||
    217258              scaledRHSError>1.0e-1||
    218               maximumDualError_>objectiveNorm_*1.0e-1) {
     259              maximumDualError_>objectiveNorm_*1.0e-1
     260              &&numberIterations_>50
     261              &&complementarityGap_>0.9*historyInfeasibility_[0]) {
    219262            handler_->message(CLP_BARRIER_EXIT2,messages_)
    220263              <<saveIteration
     
    223266          }
    224267          if (complementarityGap_>0.95*checkGap&&bestObjectiveGap<1.0e-3&&
    225               numberIterations_>saveIteration+5) {
     268              (numberIterations_>saveIteration+5||numberIterations_>100)) {
    226269            handler_->message(CLP_BARRIER_EXIT2,messages_)
    227270              <<saveIteration
     
    390433    findStepLength(phase);
    391434    nextGap=complementarityGap(nextNumber,nextNumberItems,1);
     435    debugMove(0,actualPrimalStep_,actualDualStep_);
     436    //debugMove(0,1.0e-7,1.0e-7);
    392437    double affineGap=nextGap;
    393438    int bestPhase=0;
     
    395440    // ?
    396441    bestNextGap=max(nextGap,0.8*complementarityGap_);
    397     bestNextGap=max(nextGap,0.99*complementarityGap_);
     442    //bestNextGap=max(nextGap,0.99*complementarityGap_);
    398443    if (complementarityGap_>1.0e-4*numberComplementarityPairs_) {
    399444      //std::cout <<"predicted duality gap "<<nextGap<<std::endl;
     
    413458    } else {
    414459      double phi;
    415       if (numberComplementarityPairs_<=500) {
     460      if (numberComplementarityPairs_<=5000) {
    416461        phi=pow((double) numberComplementarityPairs_,2.0);
    417462      } else {
     
    425470    //save information
    426471    double product=affineProduct();
    427     //#define ALWAYS
    428 #ifndef ALWAYS
    429472#if 0
    430473    //can we do corrector step?
     
    456499    }
    457500#endif
    458     if (complementarityGap_*(beta2-tau)+product-mu_*numberComplementarityPairs_<0.0) {
     501    if (complementarityGap_*(beta2-tau)+product-mu_*numberComplementarityPairs_<0.0&&0) {
    459502#ifdef SOME_DEBUG
    460503      printf("failed 1 product %.18g mu %.18g - %.18g < 0.0, nextGap %.18g\n",product,mu_,
     
    476519      phase=1;
    477520    }
    478 #else
    479     phase=1;
     521    // If bad gap - try standard primal dual
     522    if (nextGap>complementarityGap_*1.0001)
     523      goodMove=false;
     524    if ((modeSwitch&2)!=0)
     525      goodMove=false;
     526    if (goodMove&&doCorrector) {
     527      memcpy(saveX,deltaX_,numberTotal*sizeof(double));
     528      memcpy(saveY,deltaY_,numberRows_*sizeof(double));
     529      memcpy(saveZ,deltaZ_,numberTotal*sizeof(double));
     530      memcpy(saveT,deltaT_,numberTotal*sizeof(double));
     531      memcpy(saveSL,deltaSL_,numberTotal*sizeof(double));
     532      memcpy(saveSU,deltaSU_,numberTotal*sizeof(double));
     533#ifdef HALVE
     534      double savePrimalStep = actualPrimalStep_;
     535      double saveDualStep = actualDualStep_;
     536      double saveMu = mu_;
    480537#endif
    481     if (goodMove&&doCorrector) {
    482538      //set up for next step
    483539      setupForSolve(phase);
     
    502558               directionAccuracy2,testValue);
    503559#endif
    504         bestNextGap=COIN_DBL_MAX;
    505560      }
    506561      if (goodMove) {
     562        phase=1;
    507563        findStepLength(phase);
    508564        nextGap = complementarityGap(nextNumber,nextNumberItems,1);
    509 #ifndef ALWAYS
    510         if (numberIterations_>=-77) {
    511           goodMove=checkGoodMove(true,bestNextGap);
    512           if (!goodMove) {
     565        debugMove(1,actualPrimalStep_,actualDualStep_);
     566        //debugMove(1,1.0e-7,1.0e-7);
     567        goodMove=checkGoodMove(true,bestNextGap,allowIncreasingGap);
     568        if (!goodMove) {
    513569#ifdef SOME_DEBUG
    514             printf("checkGoodMove failed\n");
     570          printf("checkGoodMove failed\n");
    515571#endif
    516             if ((affineGap<0.5*complementarityGap_&&complementarityGap_>0.9*checkGap)
    517               ||(affineGap<0.95*complementarityGap_&&
    518                  complementarityGap_>0.99*checkGap)) {
    519               // Back to affine
    520               phase=0;
    521               // Try primal dual step instead - but with small mu
    522               if (numberIterations_<100)
    523                 phase=2;
    524               double floatNumber;
    525               floatNumber = 2.0*numberComplementarityPairs_;
    526               mu_=complementarityGap_/floatNumber;
    527               double mu1=mu_;
    528               double phi;
    529               if (numberComplementarityPairs_<=500) {
    530                 phi=pow((double) numberComplementarityPairs_,2.0);
    531               } else {
    532                 phi=pow((double) numberComplementarityPairs_,1.5);
    533                 if (phi<500.0*500.0) {
    534                   phi=500.0*500.0;
    535                 }
    536               }
    537               mu_=complementarityGap_/phi;
    538               //printf("pd mu %g, alternate %g, smallest %g\n",
    539               //     mu_,mu1,smallestPrimalDualMu);
    540               mu_ = sqrt(mu_*mu1);
    541               mu_=mu1*0.8;
    542               if (numberIterations_>100)
    543                 mu_ *=0.1; // so close to affine
    544               setupForSolve(phase);
    545               directionAccuracy=findDirectionVector(phase);
    546               findStepLength(phase);
    547               nextGap=complementarityGap(nextNumber,nextNumberItems,1);
    548               goodMove=true;
    549             }
    550           }
    551         } else {
    552           goodMove=true;
    553         }
     572        }
     573      }
     574#ifdef HALVE
     575      int nHalve=0;
     576      // relax test
     577      bestNextGap=max(bestNextGap,0.9*complementarityGap_);
     578      while (!goodMove) {
     579        mu_=saveMu;
     580        actualPrimalStep_ = savePrimalStep;
     581        actualDualStep_ = saveDualStep;
     582        int i;
     583        //printf("halve %d\n",nHalve);
     584        nHalve++;
     585        const double lambda=0.5;
     586        for (i=0;i<numberRows_;i++)
     587          deltaY_[i] = lambda*deltaY_[i]+(1.0-lambda)*saveY[i];
     588        for (i=0;i<numberTotal;i++) {
     589          deltaX_[i] = lambda*deltaX_[i]+(1.0-lambda)*saveX[i];
     590          deltaZ_[i] = lambda*deltaZ_[i]+(1.0-lambda)*saveZ[i];
     591          deltaT_[i] = lambda*deltaT_[i]+(1.0-lambda)*saveT[i];
     592          deltaSL_[i] = lambda*deltaSL_[i]+(1.0-lambda)*saveSL[i];
     593          deltaSU_[i] = lambda*deltaSU_[i]+(1.0-lambda)*saveSU[i];
     594        }
     595        //memcpy(deltaX_,saveX,numberTotal*sizeof(double));
     596        //memcpy(deltaY_,saveY,numberRows_*sizeof(double));
     597        //memcpy(deltaZ_,saveZ,numberTotal*sizeof(double));
     598        //memcpy(deltaT_,saveT,numberTotal*sizeof(double));
     599        //memcpy(deltaSL_,saveSL,numberTotal*sizeof(double));
     600        //memcpy(deltaSU_,saveSU,numberTotal*sizeof(double));
     601        findStepLength(1);
     602        nextGap = complementarityGap(nextNumber,nextNumberItems,1);
     603        goodMove=checkGoodMove(true,bestNextGap,allowIncreasingGap);
     604        if (nHalve>10)
     605          break;
     606        //assert (goodMove);
     607      }
     608      if (nHalve&&handler_->logLevel()>1)
     609        printf("halved %d times\n",nHalve);
    554610#endif
    555       }
    556611    }
    557612    //bestPhase=-1;
     
    560615      // Just primal dual step
    561616      double floatNumber;
    562       floatNumber = 2.5*numberComplementarityPairs_;
    563       //floatNumber = 1.5*numberComplementarityItems_;
     617      floatNumber = 2.0*numberComplementarityPairs_;
     618      //floatNumber = numberComplementarityItems_;
    564619      double saveMu=mu_; // use one from predictor corrector
    565620      mu_=complementarityGap_/floatNumber;
     621      // If going well try small mu
     622      mu_ *= sqrt((1.0-lastStep)/(1.0+10.0*lastStep));
    566623      double mu1=mu_;
    567624      double phi;
     
    588645      findStepLength(2);
    589646      // just for debug
     647      bestNextGap = complementarityGap_;
    590648      nextGap=complementarityGap(nextNumber,nextNumberItems,2);
    591       if (nextGap>0.9*complementarityGap_&&bestPhase==0&&affineGap<nextGap) {
     649      debugMove(2,actualPrimalStep_,actualDualStep_);
     650      //debugMove(2,1.0e-7,1.0e-7);
     651      checkGoodMove(false, bestNextGap,allowIncreasingGap);
     652      if (nextGap>0.9*complementarityGap_&&bestPhase==0&&affineGap<nextGap
     653          &&numberIterations_>80) {
    592654        // Back to affine
    593655        phase=0;
    594         // no
    595         if (numberIterations_<80)
    596           phase=2;
    597         mu_ *= 0.5;
    598656        setupForSolve(phase);
    599657        directionAccuracy=findDirectionVector(phase);
    600658        findStepLength(phase);
    601659        nextGap=complementarityGap(nextNumber,nextNumberItems,1);
     660        bestNextGap = complementarityGap_;
     661        //checkGoodMove(false, bestNextGap,allowIncreasingGap);
    602662      }
    603663    }
     
    606666    if (!goodMove)
    607667      mu_=nextGap / ((double) 1.1*nextNumber);
    608     goodMove=true;
     668    //goodMove=true; // don't do if was bad move
    609669    //goodMove=false; //TEMP
    610670    // Do centering steps
     
    616676    if (actualDualStep_>0.9&&actualPrimalStep_>0.9)
    617677      goodMove=false; // don't bother
     678    if ((modeSwitch&1)!=0)
     679      goodMove=false;
    618680    while (goodMove&&numberTries<5) {
    619681      goodMove=false;
     
    629691      findDirectionVector(3);
    630692      findStepLength(3);
     693      debugMove(3,actualPrimalStep_,actualDualStep_);
     694      //debugMove(3,1.0e-7,1.0e-7);
    631695      double xGap = complementarityGap(nextNumber,nextNumberItems,3);
    632696      // If one small then that's the one that counts
     
    638702        checkPrimal=2.0*checkDual;
    639703      }
    640       if (actualPrimalStep_<=checkPrimal||
    641           actualDualStep_<=checkDual||
     704      if (actualPrimalStep_<checkPrimal||
     705          actualDualStep_<checkDual||
    642706          (xGap>nextGap&&xGap>0.9*complementarityGap_)) {
    643707        //if (actualPrimalStep_<=checkPrimal||
     
    685749  delete [] saveZ;
    686750  delete [] saveT;
     751  delete [] saveSL;
     752  delete [] saveSU;
    687753  if (savePi) {
    688754    //std::cout<<"Restoring from iteration "<<saveIteration<<std::endl;
     
    700766  multiplyAdd(NULL,numberTotal,0.0,cost_,scaleFactor_);
    701767  multiplyAdd(NULL,numberRows_,0,dual_,scaleFactor_);
    702   CoinMemcpyN(cost_,numberColumns_,reducedCost_);
    703   matrix_->transposeTimes(-1.0,dual_,reducedCost_);
     768  checkSolution();
    704769  CoinMemcpyN(reducedCost_,numberColumns_,dj_);
    705   checkSolution();
     770  // Restore quadratic objective if necessary
     771  if (saveObjective) {
     772    delete objective_;
     773    objective_=saveObjective;
     774  }
    706775  handler_->message(CLP_BARRIER_END,messages_)
    707776    <<sumPrimalInfeasibilities_
     
    777846            chosenPrimalSequence=iColumn;
    778847          } else {
    779             printf("small %d delta %g newZ %g step %g\n",iColumn,delta,newZ,newStep);
     848            //printf("small %d delta %g newZ %g step %g\n",iColumn,delta,newZ,newStep);
    780849          }
    781850        }
     
    797866            chosenPrimalSequence=iColumn;
    798867          } else {
    799             printf("small %d delta %g newT %g step %g\n",iColumn,delta,newT,newStep);
     868            //printf("small %d delta %g newT %g step %g\n",iColumn,delta,newT,newStep);
    800869          }
    801870        }
     
    815884    actualDualStep_=1.0;
    816885  }
     886  // See if quadratic objective
     887#ifndef NO_RTTI
     888  ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
     889#else
     890  ClpQuadraticObjective * quadraticObj = NULL;
     891  if (objective_->type()==2)
     892    quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
     893#endif
     894  if (quadraticObj) {
     895    // Use smaller
     896    actualDualStep_=min(actualDualStep_,actualPrimalStep_);
     897    actualPrimalStep_=actualDualStep_;
     898  }
    817899#ifdef FULL_DEBUG
    818900  if (phase==3){
     
    822904      if (!flagged(iColumn)) {
    823905        if (lowerBound(iColumn)) {
    824           double change = rhsL_[iColumn] + deltaX_[iColumn];
     906          double change = -rhsL_[iColumn] + deltaX_[iColumn];
    825907          double dualValue=zVec[iColumn]+actualDualStep_*deltaZ_[iColumn];
    826908          double primalValue=lowerSlack[iColumn]+actualPrimalStep_*change;
     
    843925  }
    844926#endif
    845 #ifdef SOME_DEBUG
     927#ifdef SOME_DEBUG_not
    846928  {
    847929    double largestL=0.0;
     
    856938      if (!flagged(iColumn)) {
    857939        if (lowerBound(iColumn)) {
    858           double change = rhsL_[iColumn] + deltaX_[iColumn];
     940          double change = -rhsL_[iColumn] + deltaX_[iColumn];
    859941          double dualValue=zVec[iColumn]+actualDualStep_*deltaZ_[iColumn];
    860942          double primalValue=lowerSlack[iColumn]+actualPrimalStep_*change;
     
    888970      if (!flagged(iColumn)) {
    889971        if (lowerBound(iColumn)) {
    890           double change = rhsL_[iColumn] + deltaX_[iColumn];
     972          double change = -rhsL_[iColumn] + deltaX_[iColumn];
    891973          double dualValue=zVec[iColumn]+actualDualStep_*deltaZ_[iColumn];
    892974          double primalValue=lowerSlack[iColumn]+actualPrimalStep_*change;
     
    9241006  return directionNorm;
    9251007}
    926 //#define KKT 2
    9271008/* Does solve. region1 is for deltaX (columns+rows), region2 for deltaPi (rows) */
    9281009void
     
    9431024  }
    9441025  int iColumn;
    945 #if KKT<2
    946   for (iColumn=0;iColumn<numberTotal;iColumn++)
    947     region1[iColumn] = region1In[iColumn]*diagonal_[iColumn];
    948   multiplyAdd(region1+numberColumns_,numberRows_,-1.0,region2,1.0);
    949   matrix_->times(1.0,region1,region2);
    950   double maximumRHS = maximumAbsElement(region2,numberRows_);
    951   double scale=1.0;
    952   double unscale=1.0;
    953   if (maximumRHS>1.0e-30) {
    954     if (maximumRHS<=0.5) {
    955       double factor=2.0;
    956       while (maximumRHS<=0.5) {
    957         maximumRHS*=factor;
    958         scale*=factor;
    959       } /* endwhile */
    960     } else if (maximumRHS>=2.0&&maximumRHS<=COIN_DBL_MAX) {
    961       double factor=0.5;
    962       while (maximumRHS>=2.0) {
    963         maximumRHS*=factor;
    964         scale*=factor;
    965       } /* endwhile */
    966     }
    967     unscale=diagonalScaleFactor_/scale;
     1026  if (cholesky_->type()<20) {
     1027    // not KKT
     1028    for (iColumn=0;iColumn<numberTotal;iColumn++)
     1029      region1[iColumn] = region1In[iColumn]*diagonal_[iColumn];
     1030    multiplyAdd(region1+numberColumns_,numberRows_,-1.0,region2,1.0);
     1031    matrix_->times(1.0,region1,region2);
     1032    double maximumRHS = maximumAbsElement(region2,numberRows_);
     1033    double scale=1.0;
     1034    double unscale=1.0;
     1035    if (maximumRHS>1.0e-30) {
     1036      if (maximumRHS<=0.5) {
     1037        double factor=2.0;
     1038        while (maximumRHS<=0.5) {
     1039          maximumRHS*=factor;
     1040          scale*=factor;
     1041        } /* endwhile */
     1042      } else if (maximumRHS>=2.0&&maximumRHS<=COIN_DBL_MAX) {
     1043        double factor=0.5;
     1044        while (maximumRHS>=2.0) {
     1045          maximumRHS*=factor;
     1046          scale*=factor;
     1047        } /* endwhile */
     1048      }
     1049      unscale=diagonalScaleFactor_/scale;
     1050    } else {
     1051      //effectively zero
     1052      scale=0.0;
     1053      unscale=0.0;
     1054    }
     1055    multiplyAdd(NULL,numberRows_,0.0,region2,scale);
     1056    cholesky_->solve(region2);
     1057    multiplyAdd(NULL,numberRows_,0.0,region2,unscale);
     1058    multiplyAdd(region2,numberRows_,-1.0,region1+numberColumns_,0.0);
     1059    CoinZeroN(region1,numberColumns_);
     1060    matrix_->transposeTimes(1.0,region2,region1);
     1061    for (iColumn=0;iColumn<numberTotal;iColumn++)
     1062      region1[iColumn] = (region1[iColumn]-region1In[iColumn])*diagonal_[iColumn];
    9681063  } else {
    969     //effectively zero
    970     scale=0.0;
    971     unscale=0.0;
    972   }
    973   multiplyAdd(NULL,numberRows_,0.0,region2,scale);
    974   cholesky_->solve(region2);
    975   multiplyAdd(NULL,numberRows_,0.0,region2,unscale);
    976   multiplyAdd(region2,numberRows_,-1.0,region1+numberColumns_,0.0);
    977   CoinZeroN(region1,numberColumns_);
    978   matrix_->transposeTimes(1.0,region2,region1);
    979   for (iColumn=0;iColumn<numberTotal;iColumn++)
    980     region1[iColumn] = (region1[iColumn]-region1In[iColumn])*diagonal_[iColumn];
    981 #else
    982   for (iColumn=0;iColumn<numberTotal;iColumn++)
    983     region1[iColumn] = region1In[iColumn];
    984   cholesky_->solveKKT(region1,region2,diagonal_,diagonalScaleFactor_);
    985 #endif
     1064    for (iColumn=0;iColumn<numberTotal;iColumn++)
     1065      region1[iColumn] = region1In[iColumn];
     1066    cholesky_->solveKKT(region1,region2,diagonal_,diagonalScaleFactor_);
     1067  }
    9861068  if (saveRegion2) {
    9871069    //refine?
     
    10151097  //if flagged then entries zero so can do
    10161098  // For KKT separate out
    1017 #ifndef KKT
     1099  double * region1Save=NULL;//for refinement
    10181100  int iColumn;
    1019   for (iColumn=0;iColumn<numberTotal;iColumn++)
    1020     deltaX_[iColumn] = workArray[iColumn] - solution_[iColumn];
    1021   multiplyAdd(deltaX_+numberColumns_,numberRows_,-1.0,deltaY_,0.0);
    1022   matrix_->times(1.0,deltaX_,deltaY_);
    1023 #else
    1024   // regions in will be workArray and newError
    1025   // regions out deltaX_ and deltaY_
    1026   multiplyAdd(solution_+numberColumns_,numberRows_,1.0,newError,0.0);
    1027   matrix_->times(-1.0,solution_,newError);
    1028   double * region1Save=NULL;//for refinement
    1029 #endif
     1101  if (cholesky_->type()<20) {
     1102    int iColumn;
     1103    for (iColumn=0;iColumn<numberTotal;iColumn++)
     1104      deltaX_[iColumn] = workArray[iColumn] - solution_[iColumn];
     1105    multiplyAdd(deltaX_+numberColumns_,numberRows_,-1.0,deltaY_,0.0);
     1106    matrix_->times(1.0,deltaX_,deltaY_);
     1107  } else {
     1108    // regions in will be workArray and newError
     1109    // regions out deltaX_ and deltaY_
     1110    multiplyAdd(solution_+numberColumns_,numberRows_,1.0,newError,0.0);
     1111    matrix_->times(-1.0,solution_,newError);
     1112    // This is inefficient but just for now get values which will be in deltay
     1113    int iColumn;
     1114    for (iColumn=0;iColumn<numberTotal;iColumn++)
     1115      deltaX_[iColumn] = workArray[iColumn] - solution_[iColumn];
     1116    multiplyAdd(deltaX_+numberColumns_,numberRows_,-1.0,deltaY_,0.0);
     1117    matrix_->times(1.0,deltaX_,deltaY_);
     1118  }
    10301119  bool goodSolve=false;
    10311120  double * regionSave=NULL;//for refinement
     
    10361125    double lastError=relativeError;
    10371126    goodSolve=true;
    1038     double maximumRHS = maximumAbsElement(deltaY_,numberRows_);
    1039     double saveMaximum = maximumRHS;
    1040 #ifndef KKT
    1041     double scale=1.0;
    1042     double unscale=1.0;
    1043     if (maximumRHS>1.0e-30) {
    1044       if (maximumRHS<=0.5) {
    1045         double factor=2.0;
    1046         while (maximumRHS<=0.5) {
    1047           maximumRHS*=factor;
    1048           scale*=factor;
    1049         } /* endwhile */
    1050       } else if (maximumRHS>=2.0&&maximumRHS<=COIN_DBL_MAX) {
    1051         double factor=0.5;
    1052         while (maximumRHS>=2.0) {
    1053           maximumRHS*=factor;
    1054           scale*=factor;
    1055         } /* endwhile */
    1056       }
    1057       unscale=diagonalScaleFactor_/scale;
     1127    double maximumRHS;
     1128    double saveMaximum;
     1129    maximumRHS = maximumAbsElement(deltaY_,numberRows_);
     1130    saveMaximum = maximumRHS;
     1131    if (cholesky_->type()<20) {
     1132      // no kkt
     1133      double scale=1.0;
     1134      double unscale=1.0;
     1135      if (maximumRHS>1.0e-30) {
     1136        if (maximumRHS<=0.5) {
     1137          double factor=2.0;
     1138          while (maximumRHS<=0.5) {
     1139            maximumRHS*=factor;
     1140            scale*=factor;
     1141          } /* endwhile */
     1142        } else if (maximumRHS>=2.0&&maximumRHS<=COIN_DBL_MAX) {
     1143          double factor=0.5;
     1144          while (maximumRHS>=2.0) {
     1145            maximumRHS*=factor;
     1146            scale*=factor;
     1147          } /* endwhile */
     1148        }
     1149        unscale=diagonalScaleFactor_/scale;
     1150      } else {
     1151        //effectively zero
     1152        scale=0.0;
     1153        unscale=0.0;
     1154      }
     1155      //printf("--putting scales to 1.0\n");
     1156      //scale=1.0;
     1157      //unscale=1.0;
     1158      multiplyAdd(NULL,numberRows_,0.0,deltaY_,scale);
     1159      cholesky_->solve(deltaY_);
     1160      multiplyAdd(NULL,numberRows_,0.0,deltaY_,unscale);
     1161#if 0
     1162      {
     1163        printf("deltay\n");
     1164        for (int i=0;i<numberRows_;i++)
     1165          printf("%d %.18g\n",i,deltaY_[i]);
     1166      }
     1167      exit(66);
     1168#endif
     1169      if (numberTries) {
     1170        //refine?
     1171        double scaleX=1.0;
     1172        if (lastError>1.0e-5)
     1173          scaleX=0.8;
     1174        multiplyAdd(regionSave,numberRows_,1.0,deltaY_,scaleX);
     1175      }
     1176      //CoinZeroN(newError,numberRows_);
     1177      multiplyAdd(deltaY_,numberRows_,-1.0,deltaX_+numberColumns_,0.0);
     1178      CoinZeroN(deltaX_,numberColumns_);
     1179      matrix_->transposeTimes(1.0,deltaY_,deltaX_);
     1180      //if flagged then entries zero so can do
     1181      for (iColumn=0;iColumn<numberTotal;iColumn++)
     1182        deltaX_[iColumn] = deltaX_[iColumn]*diagonal_[iColumn]
     1183          -workArray[iColumn];
    10581184    } else {
    1059       //effectively zero
    1060       scale=0.0;
    1061       unscale=0.0;
    1062     }
    1063     //printf("--putting scales to 1.0\n");
    1064     //scale=1.0;
    1065     //unscale=1.0;
    1066     multiplyAdd(NULL,numberRows_,0.0,deltaY_,scale);
    1067     cholesky_->solve(deltaY_);
    1068     multiplyAdd(NULL,numberRows_,0.0,deltaY_,unscale);
    1069 #if 0
    1070  {
    1071    printf("deltay\n");
    1072    for (int i=0;i<numberRows_;i++)
    1073      printf("%d %.18g\n",i,deltaY_[i]);
    1074    }
    1075     exit(66);
    1076 #endif
    1077     if (numberTries) {
    1078       //refine?
    1079       double scaleX=1.0;
    1080       if (lastError>1.0e-5)
    1081         scaleX=0.8;
    1082       multiplyAdd(regionSave,numberRows_,1.0,deltaY_,scaleX);
    1083     }
    1084     //CoinZeroN(newError,numberRows_);
    1085     multiplyAdd(deltaY_,numberRows_,-1.0,deltaX_+numberColumns_,0.0);
    1086     CoinZeroN(deltaX_,numberColumns_);
    1087     matrix_->transposeTimes(1.0,deltaY_,deltaX_);
    1088     //if flagged then entries zero so can do
    1089     for (iColumn=0;iColumn<numberTotal;iColumn++)
    1090       deltaX_[iColumn] = deltaX_[iColumn]*diagonal_[iColumn]
    1091         -workArray[iColumn];
    1092 #else
    1093     solveSystem(deltaX_, deltaY_,
    1094                 workArray,newError,region1Save,regionSave,lastError>1.0e-5);
    1095 #endif
     1185      // KKT
     1186      solveSystem(deltaX_, deltaY_,
     1187                  workArray,newError,region1Save,regionSave,lastError>1.0e-5);
     1188    }
    10961189    multiplyAdd(deltaX_+numberColumns_,numberRows_,-1.0,newError,0.0);
    10971190    matrix_->times(1.0,deltaX_,newError);
     
    11501243        if (!regionSave) {
    11511244          regionSave = new double [numberRows_];
    1152 #ifdef KKT
    1153           region1Save = new double [numberTotal];
    1154 #endif
     1245          if (cholesky_->type()>=20)
     1246            region1Save = new double [numberTotal];
    11551247        }
    11561248        CoinMemcpyN(deltaY_,numberRows_,regionSave);
    1157 #ifndef KKT
    1158         multiplyAdd(newError,numberRows_,-1.0,deltaY_,0.0);
    1159 #else
    1160         CoinMemcpyN(deltaX_,numberTotal,region1Save);
    1161         // and back to input region
    1162         CoinMemcpyN(deltaY_,numberRows_,newError);
    1163 #endif
     1249        if (cholesky_->type()<20) {
     1250          // not KKT
     1251          multiplyAdd(newError,numberRows_,-1.0,deltaY_,0.0);
     1252        } else {
     1253          // KKT
     1254          CoinMemcpyN(deltaX_,numberTotal,region1Save);
     1255          // and back to input region
     1256          CoinMemcpyN(deltaY_,numberRows_,newError);
     1257        }
    11641258      }
    11651259    } else {
     
    11671261      //bring back previous
    11681262      relativeError=lastError;
    1169       CoinMemcpyN(regionSave,numberRows_,deltaY_);
    1170 #ifndef KKT
    1171       multiplyAdd(deltaY_,numberRows_,-1.0,deltaX_+numberColumns_,0.0);
    1172       CoinZeroN(deltaX_,numberColumns_);
    1173       matrix_->transposeTimes(1.0,deltaY_,deltaX_);
    1174       //if flagged then entries zero so can do
    1175       for (iColumn=0;iColumn<numberTotal;iColumn++)
    1176         deltaX_[iColumn] = deltaX_[iColumn]*diagonal_[iColumn]
    1177           -workArray[iColumn];
    1178 #else
    1179         CoinMemcpyN(region1Save,numberTotal,deltaX_);
    1180 #endif
    1181     }
     1263      if (regionSave) {
     1264        CoinMemcpyN(regionSave,numberRows_,deltaY_);
     1265        if (cholesky_->type()<20) {
     1266          // not KKT
     1267          multiplyAdd(deltaY_,numberRows_,-1.0,deltaX_+numberColumns_,0.0);
     1268          CoinZeroN(deltaX_,numberColumns_);
     1269          matrix_->transposeTimes(1.0,deltaY_,deltaX_);
     1270          //if flagged then entries zero so can do
     1271          for (iColumn=0;iColumn<numberTotal;iColumn++)
     1272            deltaX_[iColumn] = deltaX_[iColumn]*diagonal_[iColumn]
     1273              -workArray[iColumn];
     1274        } else {
     1275          // KKT
     1276          CoinMemcpyN(region1Save,numberTotal,deltaX_);
     1277        }
     1278      } else {
     1279        // disaster
     1280        CoinFillN(deltaX_,numberTotal,1.0);
     1281        CoinFillN(deltaY_,numberRows_,1.0);
     1282        printf("bad cholesky\n");
     1283      }
     1284    }
    11821285  } /* endwhile */
    11831286  delete [] regionSave;
    1184 #ifdef KKT
    11851287  delete [] region1Save;
    1186   int iColumn;
    1187 #endif
    11881288  delete [] newError;
    11891289  // now rest
     
    11931293  double * upperSlack = upperSlack_;
    11941294  double extra=eExtra;
     1295  //multiplyAdd(deltaY_,numberRows_,1.0,deltaT_+numberColumns_,0.0);
     1296  //CoinZeroN(deltaT_,numberColumns_);
     1297  //matrix_->transposeTimes(-1.0,deltaY_,deltaT_);
    11951298
    11961299  for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
     
    11981301    deltaSL_[iColumn]=0.0;
    11991302    deltaZ_[iColumn]=0.0;
     1303    double dd=deltaT_[iColumn];
    12001304    deltaT_[iColumn]=0.0;
    12011305    if (!flagged(iColumn)) {
    12021306      double deltaX = deltaX_[iColumn];
    12031307      if (lowerBound(iColumn)) {
    1204         double deltaSL = rhsL_[iColumn]+deltaX;
     1308        double zValue = rhsZ_[iColumn];
     1309        double gHat = zValue + zVec[iColumn]*rhsL_[iColumn];
    12051310        double slack = lowerSlack[iColumn]+extra;
    1206         deltaSL_[iColumn]=deltaSL;
    1207         deltaZ_[iColumn]=(rhsZ_[iColumn]-zVec[iColumn]*deltaSL)/slack;
     1311        deltaSL_[iColumn] = -rhsL_[iColumn]+deltaX;
     1312        deltaZ_[iColumn]=(gHat-zVec[iColumn]*deltaX)/slack;
    12081313      }
    12091314      if (upperBound(iColumn)) {
    1210         double deltaSU = rhsU_[iColumn]-deltaX;
     1315        double tValue = rhsT_[iColumn];
     1316        double hHat = tValue - tVec[iColumn]*rhsU_[iColumn];
    12111317        double slack = upperSlack[iColumn]+extra;
    1212         deltaSU_[iColumn]=deltaSU;
    1213         deltaT_[iColumn]=(rhsT_[iColumn]-tVec[iColumn]*deltaSU)/slack;
     1318        deltaSU_[iColumn] = rhsU_[iColumn]-deltaX;
     1319        deltaT_[iColumn]=(hHat+tVec[iColumn]*deltaX)/slack;
     1320      }
     1321      if (0) {
     1322        // different way of calculating
     1323        double gamma2 = gamma_*gamma_;
     1324        double dZ=0.0;
     1325        double dT=0.0;
     1326        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;
     1332        double q = rhsC_[iColumn]+gamma2 * deltaX +dd;
     1333        if (primalR_)
     1334          q += deltaX*primalR_[iColumn];
     1335        dT = (gHat+hHat -slackL*q + (tValue-zValue)*deltaX)/(slackL+slackU);
     1336        dZ = dT + q;
     1337        //printf("B %d old %g %g new %g %g\n",iColumn,deltaZ_[iColumn],
     1338        //deltaT_[iColumn],dZ,dT);
     1339        if (lowerBound(iColumn)) {
     1340          if (upperBound(iColumn)) {
     1341            //printf("B %d old %g %g new %g %g\n",iColumn,deltaZ_[iColumn],
     1342            //deltaT_[iColumn],dZ,dT);
     1343            deltaT_[iColumn]=dT;
     1344            deltaZ_[iColumn]=dZ;
     1345          } else {
     1346            // just lower
     1347            //printf("L %d old %g new %g\n",iColumn,deltaZ_[iColumn],
     1348            //dZ);
     1349          }
     1350        } else {
     1351          assert (upperBound(iColumn));
     1352          //printf("U %d old %g new %g\n",iColumn,deltaT_[iColumn],
     1353          //dT);
     1354        }
    12141355      }
    12151356    }
     
    12231364  int iColumn;
    12241365  double tolerance = primalTolerance();
    1225   for (iColumn=0;iColumn<numberTotal;iColumn++) {
    1226     if (upper_[iColumn]-lower_[iColumn]>tolerance)
     1366  // See if quadratic objective
     1367#ifndef NO_RTTI
     1368  ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
     1369#else
     1370  ClpQuadraticObjective * quadraticObj = NULL;
     1371  if (objective_->type()==2)
     1372    quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
     1373#endif
     1374  if (!quadraticObj) {
     1375    for (iColumn=0;iColumn<numberTotal;iColumn++) {
     1376      if (upper_[iColumn]-lower_[iColumn]>tolerance)
     1377        clearFixed(iColumn);
     1378      else
     1379        setFixed(iColumn);
     1380    }
     1381  } else {
     1382    // try leaving fixed
     1383    for (iColumn=0;iColumn<numberTotal;iColumn++)
    12271384      clearFixed(iColumn);
    1228     else
    1229       setFixed(iColumn);
    12301385  }
    12311386  double initialValue =1.0e-12;
    12321387
    1233   double maximumObjective=1.0e-12;
     1388  double maximumObjective=0.0;
    12341389  double objectiveNorm2=0.0;
    12351390  getNorms(cost_,numberTotal,maximumObjective,objectiveNorm2);
    1236   if (maximumObjective<1.0e-12) {
    1237     maximumObjective=1.0e-12;
     1391  if (!maximumObjective) {
     1392    maximumObjective=1.0; // objective all zero
    12381393  }
    12391394  objectiveNorm2=sqrt(objectiveNorm2)/(double) numberTotal;
     
    12511406    multiplyAdd(NULL,numberTotal,0.0,cost_,1.0/scaleFactor_);
    12521407    objectiveNorm_=maximumObjective/scaleFactor_;
    1253   }
     1408  }
     1409  // See if quadratic objective
     1410  if (quadraticObj) {
     1411    // If scaled then really scale matrix
     1412    double scaleFactor =
     1413      scaleFactor_*optimizationDirection_*objectiveScale_*
     1414      rhsScale_;
     1415    if ((scalingFlag_>0&&rowScale_)||scaleFactor != 1.0) {
     1416      CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
     1417      const int * columnQuadratic = quadratic->getIndices();
     1418      const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
     1419      const int * columnQuadraticLength = quadratic->getVectorLengths();
     1420      double * quadraticElement = quadratic->getMutableElements();
     1421      int numberColumns = quadratic->getNumCols();
     1422      double scale = 1.0/scaleFactor;
     1423      if (scalingFlag_&&rowScale_) {
     1424        for (int iColumn=0;iColumn<numberColumns;iColumn++) {
     1425          double scaleI = columnScale_[iColumn]*scale;
     1426          for (CoinBigIndex j=columnQuadraticStart[iColumn];
     1427               j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     1428            int jColumn = columnQuadratic[j];
     1429            double scaleJ = columnScale_[jColumn];
     1430            quadraticElement[j] *= scaleI*scaleJ;
     1431            objectiveNorm_ = max(objectiveNorm_,fabs(quadraticElement[j]));
     1432          }
     1433        }
     1434      } else {
     1435        // not scaled
     1436        for (int iColumn=0;iColumn<numberColumns;iColumn++) {
     1437          for (CoinBigIndex j=columnQuadraticStart[iColumn];
     1438               j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     1439            quadraticElement[j] *= scale;
     1440            objectiveNorm_ = max(objectiveNorm_,fabs(quadraticElement[j]));
     1441          }
     1442        }
     1443      }
     1444    }
     1445  }
    12541446  baseObjectiveNorm_=objectiveNorm_;
    12551447  //accumulate fixed in dj region (as spare)
     
    13491541  }
    13501542  delete [] rowsDropped;
    1351 #ifndef KKT
    1352   cholesky_->solve(errorRegion_);
    1353   //create information for solution
    1354   multiplyAdd(errorRegion_,numberRows_,-1.0,deltaX_+numberColumns_,0.0);
    1355   CoinZeroN(deltaX_,numberColumns_);
    1356   matrix_->transposeTimes(1.0,errorRegion_,deltaX_);
    1357 #else
    1358   // reverse sign on solution
    1359   multiplyAdd(NULL,numberRows_+numberColumns_,0.0,solution_,-1.0);
    1360   solveSystem(deltaX_,errorRegion_,solution_,NULL,NULL,NULL,false);
    1361 #endif
    1362   //do reduced costs
    1363   CoinZeroN(dj_+numberColumns_,numberRows_);
    1364   for ( iColumn=0;iColumn<numberColumns_;iColumn++) {
    1365     double value=cost_[iColumn];
    1366     if (lowerBound(iColumn)) {
    1367       value+=linearPerturbation_;
    1368     }
    1369     if (upperBound(iColumn)) {
    1370       value-=linearPerturbation_;
    1371     }
    1372     dj_[iColumn]=value;
    1373   }
    1374   initialValue=1.0e2;
    1375   if (rhsNorm_*1.0e-2>initialValue) {
    1376     initialValue=rhsNorm_*1.0e-2;
    1377   }
     1543  if (cholesky_->type()<20) {
     1544    // not KKT
     1545    cholesky_->solve(errorRegion_);
     1546    //create information for solution
     1547    multiplyAdd(errorRegion_,numberRows_,-1.0,deltaX_+numberColumns_,0.0);
     1548    CoinZeroN(deltaX_,numberColumns_);
     1549    matrix_->transposeTimes(1.0,errorRegion_,deltaX_);
     1550  } else {
     1551    // KKT
     1552    // reverse sign on solution
     1553    multiplyAdd(NULL,numberRows_+numberColumns_,0.0,solution_,-1.0);
     1554    solveSystem(deltaX_,errorRegion_,solution_,NULL,NULL,NULL,false);
     1555  }
     1556  //initialValue=1.0e2;
     1557  //if (rhsNorm_*1.0e-2>initialValue) {
     1558  //initialValue=rhsNorm_*1.0e-2;
     1559  //}
     1560  initialValue = max(1.0,rhsNorm_);
    13781561  double smallestBoundDifference=COIN_DBL_MAX;
    13791562  double * fakeSolution = deltaX_;
     
    13951578    <<initialValue<<objectiveNorm_
    13961579    <<CoinMessageEol;
    1397   int strategy=0;
    13981580  double extra=1.0e-10;
    13991581  double largeGap=1.0e15;
     
    14081590  double zwLarge=1.0e2*initialValue;
    14091591  //zwLarge=1.0e40;
     1592  if (cholesky_->choleskyCondition()<0.0&&cholesky_->type()<20) {
     1593    // looks bad - play safe
     1594    initialValue *=10.0;
     1595    safeObjectiveValue *= 10.0;
     1596    safeFree *= 10.0;
     1597  }
     1598  double gamma2 = gamma_*gamma_; // gamma*gamma will be added to diagonal
     1599  // First do primal side
    14101600  for ( iColumn=0;iColumn<numberTotal;iColumn++) {
    14111601    if (!flagged(iColumn)) {
    14121602      double lowerValue=lower_[iColumn];
    14131603      double upperValue=upper_[iColumn];
    1414       double objectiveValue = cost_[iColumn];
    14151604      double newValue;
    1416       double low;
    1417       double high;
    1418       double randomZ =0.1*CoinDrand48()+0.9;
    1419       double randomW =0.1*CoinDrand48()+0.9;
    14201605      double setPrimal=initialValue;
    1421       if (strategy==0) {
    1422         randomZ=1.0;
    1423         randomW=1.0;
    1424       } 
     1606      if (quadraticObj) {
     1607        // perturb primal solution a bit
     1608        //fakeSolution[iColumn]  *= 0.002*CoinDrand48()+0.999;
     1609      }
    14251610      if (lowerBound(iColumn)) {
    14261611        if (upperBound(iColumn)) {
     
    14351620            }
    14361621            newValue=fakeValue;
    1437             low=fakeValue-lowerValue;
    1438             high=upperValue-fakeValue;
    14391622          } else {
    14401623            newValue=0.5*(upperValue+lowerValue);
    1441             low=setPrimal;
    1442             high=setPrimal;
     1624          }
     1625        } else {
     1626          //just lower bound
     1627          double fakeValue=fakeSolution[iColumn];
     1628          if (fakeValue<lowerValue+setPrimal) {
     1629            fakeValue=lowerValue+setPrimal;
    14431630          }
    1444           low *= randomZ;
    1445           high *= randomW;
     1631          newValue=fakeValue;
     1632        }
     1633      } else {
     1634        if (upperBound(iColumn)) {
     1635          //just upper bound
     1636          double fakeValue=fakeSolution[iColumn];
     1637          if (fakeValue>upperValue-setPrimal) {
     1638            fakeValue=upperValue-setPrimal;
     1639          }
     1640          newValue=fakeValue;
     1641        } else {
     1642          //free
     1643          newValue=fakeSolution[iColumn];
     1644          if (newValue>=0.0) {
     1645            if (newValue<safeFree) {
     1646              newValue=safeFree;
     1647            }
     1648          } else {
     1649            if (newValue>-safeFree) {
     1650              newValue=-safeFree;
     1651            }
     1652          }
     1653        }
     1654      }
     1655      solution_[iColumn]=newValue;
     1656    } else {
     1657      // fixed
     1658      lowerSlack_[iColumn]=0.0;
     1659      upperSlack_[iColumn]=0.0;
     1660      solution_[iColumn]=lower_[iColumn];
     1661      zVec_[iColumn]=0.0;
     1662      tVec_[iColumn]=0.0;
     1663      diagonal_[iColumn]=0.0;
     1664    }
     1665  }
     1666  solutionNorm_ =  maximumAbsElement(solution_,numberTotal);
     1667  // Set bounds and do dj including quadratic
     1668  largeGap = max(1.0e7,1.02*solutionNorm_);
     1669  CoinPackedMatrix * quadratic = NULL;
     1670  const int * columnQuadratic = NULL;
     1671  const CoinBigIndex * columnQuadraticStart = NULL;
     1672  const int * columnQuadraticLength = NULL;
     1673  const double * quadraticElement = NULL;
     1674  if (quadraticObj) {
     1675    quadratic = quadraticObj->quadraticObjective();
     1676    columnQuadratic = quadratic->getIndices();
     1677    columnQuadraticStart = quadratic->getVectorStarts();
     1678    columnQuadraticLength = quadratic->getVectorLengths();
     1679    quadraticElement = quadratic->getElements();
     1680  }
     1681  double quadraticNorm=0.0;
     1682  for ( iColumn=0;iColumn<numberTotal;iColumn++) {
     1683    if (!flagged(iColumn)) {
     1684      double primalValue=solution_[iColumn];
     1685      double lowerValue=lower_[iColumn];
     1686      double upperValue=upper_[iColumn];
     1687      // Do dj
     1688      double reducedCost=cost_[iColumn];
     1689      if (lowerBound(iColumn)) {
     1690        reducedCost+=linearPerturbation_;
     1691      }
     1692      if (upperBound(iColumn)) {
     1693        reducedCost-=linearPerturbation_;
     1694      }
     1695      if (quadraticObj&&iColumn<numberColumns_) {
     1696        for (CoinBigIndex j=columnQuadraticStart[iColumn];
     1697             j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     1698          int jColumn = columnQuadratic[j];
     1699          double valueJ = solution_[jColumn];
     1700          double elementValue = quadraticElement[j];
     1701          reducedCost += valueJ*elementValue;
     1702        }
     1703        quadraticNorm = max(quadraticNorm,fabs(reducedCost));
     1704      }
     1705      dj_[iColumn]=reducedCost;
     1706      if (primalValue>lowerValue+largeGap&&primalValue<upperValue-largeGap) {
     1707        clearFixedOrFree(iColumn);
     1708        setLowerBound(iColumn);
     1709        setUpperBound(iColumn);
     1710        lowerValue=max(lowerValue,primalValue-largeGap);
     1711        upperValue=min(upperValue,primalValue+largeGap);
     1712        lower_[iColumn]=lowerValue;
     1713        upper_[iColumn]=upperValue;
     1714      }
     1715    }
     1716  }
     1717  safeObjectiveValue=max(safeObjectiveValue,quadraticNorm);
     1718  for ( iColumn=0;iColumn<numberTotal;iColumn++) {
     1719    if (!flagged(iColumn)) {
     1720      double primalValue=solution_[iColumn];
     1721      double lowerValue=lower_[iColumn];
     1722      double upperValue=upper_[iColumn];
     1723      double reducedCost=dj_[iColumn];
     1724      reducedCost = cost_[iColumn];
     1725      double low=0.0;
     1726      double high=0.0;
     1727      if (lowerBound(iColumn)) {
     1728        if (upperBound(iColumn)) {
     1729          //upper and lower bounds
     1730          if (upperValue-lowerValue>2.0*initialValue) {
     1731            low=primalValue-lowerValue;
     1732            high=upperValue-primalValue;
     1733          } else {
     1734            low=initialValue;
     1735            high=initialValue;
     1736          }
    14461737          double s = low+extra;
    14471738          double ratioZ;
     
    14661757          }
    14671758          //modify if long long way away from bound
    1468           if (objectiveValue>=0.0) {
    1469             zVec_[iColumn]=objectiveValue + safeObjectiveValue*ratioZ;
     1759          if (reducedCost>=0.0) {
     1760            zVec_[iColumn]=reducedCost + safeObjectiveValue*ratioZ;
    14701761            tVec_[iColumn]=safeObjectiveValue*ratioW;
    14711762          } else {
    14721763            zVec_[iColumn]=safeObjectiveValue*ratioZ;
    1473             tVec_[iColumn]=-objectiveValue + safeObjectiveValue*ratioW;
    1474           }
    1475           diagonal_[iColumn] = (t*s)/(s*tVec_[iColumn]+t*zVec_[iColumn]);
     1764            tVec_[iColumn]=-reducedCost + safeObjectiveValue*ratioW;
     1765          }
     1766          double gammaTerm = gamma2;
     1767          if (primalR_)
     1768            gammaTerm += primalR_[iColumn];
     1769          diagonal_[iColumn] = (t*s)/
     1770            (s*tVec_[iColumn]+t*zVec_[iColumn]+gammaTerm*t*s);
    14761771        } else {
    14771772          //just lower bound
    1478           double fakeValue=fakeSolution[iColumn];
    1479           if (fakeValue<lowerValue+setPrimal) {
    1480             fakeValue=lowerValue+setPrimal;
    1481           }
    1482           newValue=fakeValue;
    1483           low=fakeValue-lowerValue;
    1484           low *= randomZ;
     1773          low=primalValue-lowerValue;
    14851774          high=0.0;
    14861775          double s = low+extra;
     
    14951784            s=largeGap;
    14961785          }
    1497           if (objectiveValue>=0.0) {
    1498             zVec_[iColumn]=objectiveValue + safeObjectiveValue*ratioZ;
     1786          if (reducedCost>=0.0) {
     1787            zVec_[iColumn]=reducedCost + safeObjectiveValue*ratioZ;
    14991788            tVec_[iColumn]=0.0;
    15001789          } else {
     
    15021791            tVec_[iColumn]=0.0;
    15031792          }
    1504           diagonal_[iColumn] = s/zVec_[iColumn];
     1793          double gammaTerm = gamma2;
     1794          if (primalR_)
     1795            gammaTerm += primalR_[iColumn];
     1796          diagonal_[iColumn] = s/(zVec_[iColumn]+s*gammaTerm);
    15051797        }
    15061798      } else {
    15071799        if (upperBound(iColumn)) {
    15081800          //just upper bound
    1509           double fakeValue=fakeSolution[iColumn];
    1510           if (fakeValue>upperValue-setPrimal) {
    1511             fakeValue=upperValue-setPrimal;
    1512           }
    1513           newValue=fakeValue;
    15141801          low=0.0;
    1515           high=upperValue-fakeValue;
    1516           high*=randomW;
     1802          high=upperValue-primalValue;
    15171803          double t = high+extra;
    15181804          double ratioW;
     
    15261812            t=largeGap;
    15271813          }
    1528           if (objectiveValue>=0.0) {
     1814          if (reducedCost>=0.0) {
    15291815            zVec_[iColumn]=0.0;
    15301816            tVec_[iColumn]=safeObjectiveValue*ratioW;
    15311817          } else {
    15321818            zVec_[iColumn]=0.0;
    1533             tVec_[iColumn]=-objectiveValue + safeObjectiveValue*ratioW;
     1819            tVec_[iColumn]=-reducedCost + safeObjectiveValue*ratioW;
    15341820          }
    1535           diagonal_[iColumn] =  t/tVec_[iColumn];
    1536         } else {
    1537           //free
    1538           zVec_[iColumn]=safeObjectiveValue;
    1539           tVec_[iColumn]=safeObjectiveValue;
    1540           newValue=fakeSolution[iColumn];
    1541           if (newValue>=0.0) {
    1542             if (newValue<safeFree) {
    1543               newValue=safeFree;
    1544             }
    1545           } else {
    1546             if (newValue>-safeFree) {
    1547               newValue=-safeFree;
    1548             }
    1549           }
    1550           if (fabs(newValue)>1.0) {
    1551             diagonal_[iColumn]=fabs(newValue)/(zVec_[iColumn]+tVec_[iColumn]);
    1552           } else {
    1553             diagonal_[iColumn]=1.0/(zVec_[iColumn]+tVec_[iColumn]);
    1554           }
    1555           low=0.0;
    1556           high=0.0;
     1821          double gammaTerm = gamma2;
     1822          if (primalR_)
     1823            gammaTerm += primalR_[iColumn];
     1824          diagonal_[iColumn] =  t/(tVec_[iColumn]+t*gammaTerm);
    15571825        }
    15581826      }
    15591827      lowerSlack_[iColumn]=low;
    15601828      upperSlack_[iColumn]=high;
    1561       solution_[iColumn]=newValue;
    1562     } else {
    1563       lowerSlack_[iColumn]=0.0;
    1564       upperSlack_[iColumn]=0.0;
    1565       solution_[iColumn]=lower_[iColumn];
    1566       zVec_[iColumn]=0.0;
    1567       tVec_[iColumn]=0.0;
    1568       diagonal_[iColumn]=0.0;
    1569     }
    1570   }
    1571   solutionNorm_ =  maximumAbsElement(solution_,numberTotal);
    1572   // Set bounds
    1573   largeGap = max(1.0e7,1.02*solutionNorm_);
    1574   for ( iColumn=0;iColumn<numberTotal;iColumn++) {
    1575     if (!flagged(iColumn)) {
    1576       double newValue=solution_[iColumn];
    1577       double lowerValue=lower_[iColumn];
    1578       double upperValue=upper_[iColumn];
    1579       if (newValue>lowerValue+largeGap&&newValue<upperValue-largeGap) {
    1580         clearFixedOrFree(iColumn);
    1581         setLowerBound(iColumn);
    1582         setUpperBound(iColumn);
    1583         double objectiveValue = cost_[iColumn];
    1584         lowerValue=max(lowerValue,newValue-largeGap);
    1585         upperValue=min(upperValue,newValue+largeGap);
    1586         lower_[iColumn]=lowerValue;
    1587         upper_[iColumn]=upperValue;
    1588         double low;
    1589         double high;
    1590         double setPrimal=initialValue;
    1591         //upper and lower bounds
    1592         if (upperValue-lowerValue>2.0*setPrimal) {
    1593           low=newValue-lowerValue;
    1594           high=upperValue-newValue;
    1595         } else {
    1596           newValue=0.5*(upperValue+lowerValue);
    1597           low=setPrimal;
    1598           high=setPrimal;
    1599         }
    1600         double s = low+extra;
    1601         double ratioZ;
    1602         if (s<zwLarge) {
    1603           ratioZ=1.0;
    1604         } else {
    1605           ratioZ=sqrt(zwLarge/s);
    1606         }
    1607         double t = high+extra;
    1608         double ratioW;
    1609         if (t<zwLarge) {
    1610           ratioW=1.0;
    1611         } else {
    1612           ratioW=sqrt(zwLarge/t);
    1613         }
    1614         //modify s and t
    1615         if (s>largeGap) {
    1616           s=largeGap;
    1617         }
    1618         if (t>largeGap) {
    1619           t=largeGap;
    1620         }
    1621         //modify if long long way away from bound
    1622         if (objectiveValue>=0.0) {
    1623           zVec_[iColumn]=objectiveValue + safeObjectiveValue*ratioZ;
    1624           tVec_[iColumn]=safeObjectiveValue*ratioW;
    1625         } else {
    1626           zVec_[iColumn]=safeObjectiveValue*ratioZ;
    1627           tVec_[iColumn]=-objectiveValue + safeObjectiveValue*ratioW;
    1628         }
    1629         diagonal_[iColumn] = (t*s)/(s*tVec_[iColumn]+t*zVec_[iColumn]);
    1630         lowerSlack_[iColumn]=low;
    1631         upperSlack_[iColumn]=high;
    1632       }
    16331829    }
    16341830  }
     
    16611857  numberComplementarityPairs=0;
    16621858  numberComplementarityItems=0;
     1859  int numberTotal = numberRows_+numberColumns_;
    16631860  double toleranceGap=0.0;
    16641861  double largestGap=0.0;
     
    16851882  double * deltaT = deltaT_;
    16861883  double * deltaX = deltaX_;
    1687   for (int iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
     1884  for (int iColumn=0;iColumn<numberTotal;iColumn++) {
    16881885    if (!fixedOrFree(iColumn)) {
    16891886      numberComplementarityPairs++;
    16901887      //can collapse as if no lower bound both zVec and deltaZ 0.0
     1888      double newZ=0.0;
     1889      double newT=0.0;
    16911890      if (lowerBound(iColumn)) {
    16921891        numberComplementarityItems++;
     
    17001899          change =primal[iColumn]+deltaX[iColumn]-lowerSlack[iColumn]-lower[iColumn];
    17011900          dualValue=zVec[iColumn]+actualDualStep_*deltaZ[iColumn];
     1901          newZ=dualValue;
    17021902          primalValue=lowerSlack[iColumn]+actualPrimalStep_*change;
    17031903        }
     
    17341934          change =upper[iColumn]-primal[iColumn]-deltaX[iColumn]-upperSlack[iColumn];
    17351935          dualValue=tVec[iColumn]+actualDualStep_*deltaT[iColumn];
     1936          newT=dualValue;
    17361937          primalValue=upperSlack[iColumn]+actualPrimalStep_*change;
    17371938        }
     
    17571958      }
    17581959    }
    1759   } 
     1960  }
    17601961  if (!phase&&numberNegativeGaps) {
    17611962      handler_->message(CLP_BARRIER_NEGATIVE_GAPS,messages_)
     
    17911992  printf("phase %d in setupForSolve, mu %.18g\n",phase,mu_);
    17921993#endif
     1994  double gamma2 = gamma_*gamma_; // gamma*gamma will be added to diagonal
    17931995  switch (phase) {
    17941996  case 0:
    17951997    CoinMemcpyN(errorRegion_,numberRows_,rhsB_);
     1998    if (delta_||dualR_) {
     1999      // add in regularization
     2000      double delta2 = delta_*delta_;
     2001      for (int iRow=0;iRow<numberRows_;iRow++) {
     2002        rhsB_[iRow] -= delta2*dual_[iRow];
     2003        if (dualR_)
     2004          rhsB_[iRow] -= dualR_[iRow]*dual_[iRow];
     2005      }
     2006    }
    17962007    for (iColumn=0;iColumn<numberTotal;iColumn++) {
    17972008      rhsC_[iColumn]=0.0;
     
    18012012      rhsT_[iColumn]=0.0;
    18022013      if (!flagged(iColumn)) {
    1803         rhsC_[iColumn] = dj[iColumn]-zVec[iColumn]+tVec[iColumn];
    1804         if (lowerBound(iColumn)) {
     2014        rhsC_[iColumn] = dj[iColumn]-zVec[iColumn]+tVec[iColumn];
     2015        rhsC_[iColumn] += gamma2*solution_[iColumn];
     2016        if (primalR_)
     2017          rhsC_[iColumn] += primalR_[iColumn]*solution_[iColumn];
     2018        if (lowerBound(iColumn)) {
    18052019          rhsZ_[iColumn] = -zVec[iColumn]*(lowerSlack[iColumn]+extra);
    1806           rhsL_[iColumn] = min(0.0,primal[iColumn]-lowerSlack[iColumn]-lower[iColumn]);
    1807         }
    1808         if (upperBound(iColumn)) {
     2020          rhsL_[iColumn] = max(0.0,(lower[iColumn]+lowerSlack[iColumn])-solution_[iColumn]);
     2021        }
     2022        if (upperBound(iColumn)) {
    18092023          rhsT_[iColumn] = -tVec[iColumn]*(upperSlack[iColumn]+extra);
    1810           rhsU_[iColumn] = min(0.0,-(primal[iColumn] + upperSlack[iColumn] - upper[iColumn]));
    1811         }
     2024          rhsU_[iColumn] = min(0.0,(upper[iColumn]-upperSlack[iColumn])-primal[iColumn]);
     2025        }
    18122026      }
    18132027    }
     
    18592073    // could be stored in delta2?
    18602074    for (iColumn=0;iColumn<numberTotal;iColumn++) {
    1861       rhsC_[iColumn]=0.0;
    1862       //rhsU_[iColumn]=0.0;
    1863       //rhsL_[iColumn]=0.0;
    18642075      rhsZ_[iColumn]=0.0;
    18652076      rhsT_[iColumn]=0.0;
    18662077      if (!flagged(iColumn)) {
    1867         rhsC_[iColumn] = dj[iColumn]-zVec[iColumn]+tVec[iColumn];
    18682078        if (lowerBound(iColumn)) {
    18692079          rhsZ_[iColumn] = mu_ -zVec[iColumn]*(lowerSlack[iColumn]+extra)
    18702080            - deltaZ_[iColumn]*deltaX_[iColumn];
    18712081          // To bring in line with OSL
    1872           rhsZ_[iColumn] -= deltaZ_[iColumn]*rhsL_[iColumn];
    1873           //rhsZ_[iColumn] = mu_ -deltaZ_[iColumn]*rhsL_[iColumn]
    1874           //- deltaZ_[iColumn]*deltaX_[iColumn];
     2082          rhsZ_[iColumn] += deltaZ_[iColumn]*rhsL_[iColumn];
    18752083        }
    18762084        if (upperBound(iColumn)) {
     
    18792087          // To bring in line with OSL
    18802088          rhsT_[iColumn] -= deltaT_[iColumn]*rhsU_[iColumn];
    1881           //rhsT_[iColumn] = mu_ -deltaT_[iColumn]*rhsU_[iColumn]
    1882           //+deltaT_[iColumn]*deltaX_[iColumn];
    18832089        }
    18842090      }
     
    19222128    CoinMemcpyN(errorRegion_,numberRows_,rhsB_);
    19232129    for (iColumn=0;iColumn<numberTotal;iColumn++) {
    1924       rhsC_[iColumn]=0.0;
    19252130      rhsZ_[iColumn]=0.0;
    19262131      rhsT_[iColumn]=0.0;
    19272132      if (!flagged(iColumn)) {
    1928         rhsC_[iColumn] = dj[iColumn]-zVec[iColumn]+tVec[iColumn];
    19292133        if (lowerBound(iColumn)) {
    19302134          rhsZ_[iColumn] = mu_ - zVec[iColumn]*(lowerSlack[iColumn]+extra);
     
    19502154        if (!flagged(iColumn)) {
    19512155          if (lowerBound(iColumn)) {
    1952             double change = rhsL_[iColumn] + deltaX_[iColumn];
     2156            double change = -rhsL_[iColumn] + deltaX_[iColumn];
    19532157            double dualValue=zVec[iColumn]+dualStep*deltaZ_[iColumn];
    19542158            double primalValue=lowerSlack[iColumn]+primalStep*change;
     
    19722176              assert (value<0.0);
    19732177            }
    1974             //#define AGAIN
    1975 #ifdef AGAIN
    1976             //rhsZ_[iColumn] = mu_ -zVec[iColumn]*(lowerSlack[iColumn]+extra)
    1977             //- deltaZ_[iColumn]*deltaX_[iColumn];
    1978             // To bring in line with OSL
    1979             rhsZ_[iColumn] += deltaZ_[iColumn]*rhsL_[iColumn];
    1980 #endif
    19812178            rhsZ_[iColumn] += value;
    19822179          } 
     
    20022199              assert (value<0.0);
    20032200            }
    2004 #ifdef AGAIN
    2005             //rhsT_[iColumn] = mu_ -tVec[iColumn]*(upperSlack[iColumn]+extra)
    2006             //+deltaT_[iColumn]*deltaX_[iColumn];
    2007             // To bring in line with OSL
    2008             rhsT_[iColumn] -= deltaT_[iColumn]*rhsU_[iColumn];
    2009 #endif
    20102201            rhsT_[iColumn] += value;
    20112202          }
     
    20152206    break;
    20162207  } /* endswitch */
    2017 #ifndef KKT
    2018   for (iColumn=0;iColumn<numberTotal;iColumn++) {
    2019     double value = rhsC_[iColumn];
    2020     double zValue = rhsZ_[iColumn];
    2021     double tValue = rhsT_[iColumn];
     2208  if (cholesky_->type()<20) {
     2209    for (iColumn=0;iColumn<numberTotal;iColumn++) {
     2210      double value = rhsC_[iColumn];
     2211      double zValue = rhsZ_[iColumn];
     2212      double tValue = rhsT_[iColumn];
     2213#if 0
    20222214#if 1
    2023     if (phase==0) {
    2024       // more accurate
    2025       value = dj[iColumn];
    2026       zValue=0.0;
    2027       tValue=0.0;
    2028     } else if (phase==2) {
    2029       // more accurate
    2030       value = dj[iColumn];
    2031       zValue=mu_;
    2032       tValue=mu_;
    2033     }
     2215      if (phase==0) {
     2216        // more accurate
     2217        value = dj[iColumn];
     2218        zValue=0.0;
     2219        tValue=0.0;
     2220      } else if (phase==2) {
     2221        // more accurate
     2222        value = dj[iColumn];
     2223        zValue=mu_;
     2224        tValue=mu_;
     2225      }
    20342226#endif
    2035     assert (rhsL_[iColumn]<=0.0);
    2036     assert (rhsU_[iColumn]<=0.0);
    2037     if (lowerBound(iColumn)) {
    2038       value += (zVec[iColumn]*rhsL_[iColumn]-zValue)/
    2039         (lowerSlack[iColumn]+extra);
    2040     }
    2041     if (upperBound(iColumn)) {
    2042       value += (tValue-tVec[iColumn]*rhsU_[iColumn])/
    2043         (upperSlack[iColumn]+extra);
    2044     }
    2045     workArray[iColumn]=diagonal_[iColumn]*value;
    2046   }
     2227      assert (rhsL_[iColumn]>=0.0);
     2228      assert (rhsU_[iColumn]<=0.0);
     2229      if (lowerBound(iColumn)) {
     2230        value += (-zVec[iColumn]*rhsL_[iColumn]-zValue)/
     2231          (lowerSlack[iColumn]+extra);
     2232      }
     2233      if (upperBound(iColumn)) {
     2234        value += (tValue-tVec[iColumn]*rhsU_[iColumn])/
     2235          (upperSlack[iColumn]+extra);
     2236      }
     2237#else
     2238      if (lowerBound(iColumn)) {
     2239        double gHat = zValue + zVec[iColumn]*rhsL_[iColumn];
     2240        value -= gHat/(lowerSlack[iColumn]+extra);
     2241      }
     2242      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;
     2248    }
    20472249#if 0
    20482250    if (solution_[0]>0.0) {
     
    20552257    exit(66);
    20562258#endif
    2057 #else
    2058   for (iColumn=0;iColumn<numberTotal;iColumn++) {
    2059     double value = rhsC_[iColumn];
    2060     if (lowerBound(iColumn)) {
    2061       value -= rhsZ_[iColumn]/(lowerSlack[iColumn]+extra);
    2062       value += zVec[iColumn]*rhsL_[iColumn]/(lowerSlack[iColumn]+extra);
     2259  } else {
     2260    // KKT
     2261    for (iColumn=0;iColumn<numberTotal;iColumn++) {
     2262      double value = rhsC_[iColumn];
     2263      double zValue = rhsZ_[iColumn];
     2264      double tValue = rhsT_[iColumn];
     2265      if (lowerBound(iColumn)) {
     2266        double gHat = zValue + zVec[iColumn]*rhsL_[iColumn];
     2267        value -= gHat/(lowerSlack[iColumn]+extra);
     2268      }
     2269      if (upperBound(iColumn)) {
     2270        double hHat = tValue - tVec[iColumn]*rhsU_[iColumn];
     2271        value += hHat/(upperSlack[iColumn]+extra);
     2272      }
     2273      workArray[iColumn]=value;
    20632274    }
    2064     if (upperBound(iColumn)) {
    2065       value += rhsT_[iColumn]/(upperSlack[iColumn]+extra);
    2066       value -= tVec[iColumn]*rhsU_[iColumn]/(upperSlack[iColumn]+extra);
    2067     }
    2068       workArray[iColumn]=value;
    2069   }
    2070 #endif
     2275  }
    20712276}
    20722277//method: sees if looks plausible change in complementarity
    2073 bool ClpPredictorCorrector::checkGoodMove(const bool doCorrector,double & bestNextGap)
     2278bool ClpPredictorCorrector::checkGoodMove(const bool doCorrector,
     2279                                          double & bestNextGap,
     2280                                          bool allowIncreasingGap)
    20742281{
    20752282  const double beta3 = 0.99997;
     
    20802287  double returnGap=bestNextGap;
    20812288  double nextGap=complementarityGap(nextNumber,nextNumberItems,2);
    2082   if (nextGap>bestNextGap&&nextGap>0.9*complementarityGap_) {
     2289#ifndef NO_RTTI
     2290  ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
     2291#else
     2292  ClpQuadraticObjective * quadraticObj = NULL;
     2293  if (objective_->type()==2)
     2294    quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
     2295#endif
     2296  if (nextGap>bestNextGap&&nextGap>0.9*complementarityGap_&&doCorrector
     2297      &&!quadraticObj&&!allowIncreasingGap) {
    20832298#ifdef SOME_DEBUG
    20842299    printf("checkGood phase 1 next gap %.18g, phase 0 %.18g, old gap %.18g\n",
     
    21082323    //}
    21092324    double gap = bestNextGap;
    2110     goodMove=checkGoodMove2(step,gap);
     2325    goodMove=checkGoodMove2(step,gap,allowIncreasingGap);
    21112326    if (goodMove)
    21122327      returnGap=gap;
     
    21142329    goodMove=true;
    21152330  }
     2331  if (goodMove)
     2332    goodMove=checkGoodMove2(step,bestNextGap,allowIncreasingGap);
    21162333  if (!goodMove) {
    21172334    //try smaller of two
     
    21262343    actualPrimalStep_=step;
    21272344    actualDualStep_=step;
     2345    goodMove=checkGoodMove2(step,bestNextGap,allowIncreasingGap);
    21282346    while (!goodMove) {
    21292347      double gap = bestNextGap;
    2130       goodMove=checkGoodMove2(step,gap);
    2131       if (goodMove)
     2348      goodMove=checkGoodMove2(step,gap,allowIncreasingGap);
     2349      if (goodMove) {
    21322350        returnGap=gap;
     2351        break;
     2352      }
    21332353      if (step<1.0e-10) {
    21342354        break;
     
    22262446}
    22272447//:  checks for one step size
    2228 bool ClpPredictorCorrector::checkGoodMove2(const double move,double & bestNextGap)
     2448bool ClpPredictorCorrector::checkGoodMove2(const double move,
     2449                                           double & bestNextGap,
     2450                                           bool allowIncreasingGap)
    22292451{
    22302452  double complementarityMultiplier =1.0/numberComplementarityPairs_;
     
    22352457  int nextNumberItems;
    22362458  double nextGap=complementarityGap(nextNumber,nextNumberItems,2);
    2237   if (nextGap>bestNextGap)
     2459  if (nextGap>bestNextGap&&!allowIncreasingGap)
    22382460    return false;
    22392461  double lowerBoundGap = gamma*nextGap*complementarityMultiplier;
     
    22582480      }
    22592481      if (upperBound(iColumn)) {
    2260         double part1=upperSlack[iColumn]+actualPrimalStep_*deltaT[iColumn];
    2261         double part2=tVec[iColumn]+actualDualStep_*deltaSU[iColumn];
     2482        double part1=upperSlack[iColumn]+actualPrimalStep_*deltaSU[iColumn];
     2483        double part2=tVec[iColumn]+actualDualStep_*deltaT[iColumn];
    22622484        if (part1*part2<lowerBoundGap) {
    22632485          goodMove=false;
     
    22672489    }
    22682490  }
    2269   //      Satisfy g_p(alpha)?
     2491   double * nextDj=NULL;
     2492#ifndef NO_RTTI
     2493  ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
     2494#else
     2495  ClpQuadraticObjective * quadraticObj = NULL;
     2496  if (objective_->type()==2)
     2497    quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
     2498#endif
     2499  if (quadraticObj) {
     2500    double gamma2 = gamma_*gamma_;
     2501    nextDj = new double [numberColumns_];
     2502    double * nextSolution = new double [numberColumns_];
     2503    // put next primal into nextSolution
     2504    for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
     2505      if (!flagged(iColumn)) {
     2506        nextSolution[iColumn]=solution_[iColumn]+
     2507          actualPrimalStep_*deltaX_[iColumn];
     2508      }
     2509    }
     2510    // do reduced costs
     2511    CoinMemcpyN(cost_,numberColumns_,nextDj);
     2512    matrix_->transposeTimes(-1.0,dual_,nextDj);
     2513    matrix_->transposeTimes(-actualDualStep_,deltaY_,nextDj);
     2514    quadraticDjs(nextDj,nextSolution,1.0);
     2515    delete [] nextSolution;
     2516    CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
     2517    const int * columnQuadraticLength = quadratic->getVectorLengths();
     2518    for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
     2519      if (!fixedOrFree(iColumn)) {
     2520        double newZ=0.0;
     2521        double newT=0.0;
     2522        if (lowerBound(iColumn)) {
     2523          newZ=zVec[iColumn]+actualDualStep_*deltaZ[iColumn];
     2524        }
     2525        if (upperBound(iColumn)) {
     2526          newT=tVec[iColumn]+actualDualStep_*deltaT[iColumn];
     2527        }
     2528        if (columnQuadraticLength[iColumn]) {
     2529          double gammaTerm = gamma2;
     2530          if (primalR_)
     2531            gammaTerm += primalR_[iColumn];
     2532          double dualInfeasibility=
     2533            dj_[iColumn]-zVec[iColumn]+tVec[iColumn]
     2534            +gammaTerm*solution_[iColumn];
     2535          double newInfeasibility=
     2536            nextDj[iColumn]-newZ+newT
     2537            +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          }
     2545        }
     2546      }
     2547    }
     2548    delete [] nextDj;
     2549  }
     2550 //      Satisfy g_p(alpha)?
    22702551  if (rhsNorm_>solutionNorm_) {
    22712552    solutionNorm_=rhsNorm_;
     
    23232604    caution=true;
    23242605  }
    2325   // do reduced costs
    2326   CoinMemcpyN(dual_,numberRows_,dj_+numberColumns_);
    2327   CoinMemcpyN(cost_,numberColumns_,dj_);
    2328   matrix_->transposeTimes(-1.0,dual_,dj_);
    23292606  double extra=eExtra;
    23302607  const double largeFactor=1.0e2;
     
    23472624    qDiagonal=1.0e-8*mu_;
    23482625  }
     2626  //double nextMu = nextGap/((double)(2*numberComplementarityPairs_));
     2627  //printf("using gap of %g\n",nextMu);
    23492628  //qDiagonal *= 1.0e2;
    23502629  //largest allowable ratio of lowerSlack/zVec (etc)
     
    23842663  // When to start looking at killing (factor0
    23852664  double killFactor;
    2386   if (numberIterations_<50)
     2665  if (numberIterations_<50) {
    23872666    killFactor = 1.0;
    2388   else
     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;
     2673  } else {
    23892674    killFactor = 1.0e5;
    2390   for (int iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
     2675    stepLength_=0.999995;
     2676  }
     2677  // put next primal into deltaSL_
     2678  int iColumn;
     2679  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    }
     2685  }
     2686  // do reduced costs
     2687  CoinMemcpyN(dual_,numberRows_,dj_+numberColumns_);
     2688  CoinMemcpyN(cost_,numberColumns_,dj_);
     2689  matrix_->transposeTimes(-1.0,dual_,dj_);
     2690  double quadraticOffset=quadraticDjs(dj_,deltaSL_,1.0);
     2691  double gamma2 = gamma_*gamma_; // gamma*gamma will be added to diagonal
     2692  double gammaOffset=0.0;
     2693  for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
    23912694    if (!flagged(iColumn)) {
    23922695      double reducedCost=dual[iColumn];
     
    23992702      double oldPrimal = primal[iColumn];
    24002703      double newPrimal=primal[iColumn]+actualPrimalStep_*thisWeight;
    2401       double primalObjectiveThis=newPrimal*cost[iColumn];
    24022704      double dualObjectiveThis=0.0;
    2403       double t=extra;
    2404       double s=extra;
     2705      double sUpper=extra;
     2706      double sLower=extra;
    24052707      double kill;
    24062708      if (fabs(newPrimal)>1.0e4) {
    24072709        kill=killTolerance*1.0e-4*newPrimal;
    2408         //kill=killTolerance;
    2409         //kill*=1.0e-3;//???????
    24102710      } else {
    24112711        kill=killTolerance;
    2412         //kill*=1.0e-3;//???????
    24132712      }
    24142713      kill*=1.0e-3;//be conservative
    24152714      double smallerSlack=COIN_DBL_MAX;
    2416       double largerzw;
    2417       if (zValue>tValue) {
    2418         largerzw=zValue;
    2419       } else {
    2420         largerzw=tValue;
    2421       }
    24222715      bool fakeOldBounds=false;
    24232716      bool fakeNewBounds=false;
     
    24592752        fakeOldBounds=true;
    24602753      }
     2754      double lowerBoundInfeasibility=0.0;
     2755      double upperBoundInfeasibility=0.0;
    24612756      if (lowerBound(iColumn)) {
    24622757        double oldSlack = lowerSlack[iColumn];
     
    24722767          epsilon=1.0e-5;
    24732768        }
    2474         //for changing slack
    2475         double zValue2 = zValue;
    2476         if (zValue2<epsilonBase) {
    2477           zValue2=epsilonBase;
    2478         }
     2769        //epsilon=1.0e-14;
    24792770        //make sure reasonable
    24802771        if (zValue<epsilon) {
    24812772          zValue=epsilon;
    24822773        }
    2483         //store modified zVec
    2484         //no zVec[iColumn]=zValue;
    24852774        double feasibleSlack=newPrimal-lower[iColumn];
    24862775        if (feasibleSlack>0.0&&newSlack>0.0) {
     
    25062795          smallerSlack=newSlack;
    25072796        }
    2508         double infeasibility = fabs(newPrimal-lowerSlack[iColumn]-lower[iColumn]);
    2509         if (infeasibility>maximumBoundInfeasibility) {
    2510           maximumBoundInfeasibility=infeasibility;
    2511         }
     2797        lowerBoundInfeasibility = fabs(newPrimal-lowerSlack[iColumn]-lower[iColumn]);
    25122798        if (lowerSlack[iColumn]<=kill*killFactor&&fabs(newPrimal-lower[iColumn])<=kill*killFactor) {
    25132799          double step = min(actualPrimalStep_*1.1,1.0);
     
    25262812          //cout<<j<<" l "<<reducedCost<<" "<<zVec[iColumn]<<endl;
    25272813        } else {
    2528           s+=lowerSlack[iColumn];
     2814          sLower+=lowerSlack[iColumn];
    25292815        }
    25302816      }
    25312817      if (upperBound(iColumn)) {
    2532         //reducedCost-=perturbation;
    2533         //primalObjectiveThis-=perturbation*newPrimal;
    25342818        double oldSlack = upperSlack[iColumn];
    25352819        double newSlack;
     
    25442828          epsilon=1.0e-5;
    25452829        }
    2546         //for changing slack
    2547         double tValue2 = tValue;
    2548         if (tValue2<epsilonBase) {
    2549           tValue2=epsilonBase;
    2550         }
    25512830        //make sure reasonable
     2831        //epsilon=1.0e-14;
    25522832        if (tValue<epsilon) {
    25532833          tValue=epsilon;
    25542834        }
    2555         //store modified tVec
    2556         //no tVec[iColumn]=tValue;
    25572835        double feasibleSlack=upper[iColumn]-newPrimal;
    25582836        if (feasibleSlack>0.0&&newSlack>0.0) {
     
    25782856          smallerSlack=newSlack;
    25792857        }
    2580         double infeasibility = fabs(newPrimal+upperSlack[iColumn]-upper[iColumn]);
    2581         if (infeasibility>maximumBoundInfeasibility) {
    2582           maximumBoundInfeasibility=infeasibility;
    2583         }
     2858        upperBoundInfeasibility = fabs(newPrimal+upperSlack[iColumn]-upper[iColumn]);
    25842859        if (upperSlack[iColumn]<=kill*killFactor&&fabs(newPrimal-upper[iColumn])<=kill*killFactor) {
    25852860          double step = min(actualPrimalStep_*1.1,1.0);
     
    25972872          thisKilled=true;
    25982873        } else {
    2599           t+=upperSlack[iColumn];
     2874          sUpper+=upperSlack[iColumn];
    26002875        }
    26012876      }
     
    26052880      }
    26062881      if (!thisKilled) {
    2607         primalObjectiveValue+=primalObjectiveThis;
     2882        double gammaTerm = gamma2;
     2883        if (primalR_) {
     2884          gammaTerm += primalR_[iColumn];
     2885          quadraticOffset += newPrimal*newPrimal*primalR_[iColumn];
     2886        }
     2887        double dualInfeasibility=
     2888          reducedCost-zVec[iColumn]+tVec[iColumn]+gammaTerm*newPrimal;
     2889        if (fabs(dualInfeasibility)>dualTolerance) {
     2890#if 0
     2891          if (dualInfeasibility>0.0) {
     2892            // To improve we could reduce t and/or increase z
     2893            if (lowerBound(iColumn)) {
     2894              double complementarity =zVec_[iColumn]*lowerSlack[iColumn];
     2895              if (complementarity<nextMu) {
     2896                double change=
     2897                  min(dualInfeasibility,
     2898                      (nextMu-complementarity)/lowerSlack[iColumn]);
     2899                dualInfeasibility -= change;
     2900                printf("%d lb locomp %g - dual inf from %g to %g\n",
     2901                       iColumn,complementarity,dualInfeasibility+change,
     2902                       dualInfeasibility);
     2903                zVec_[iColumn] += change;
     2904                zValue = max(zVec_[iColumn],1.0e-12);
     2905              }
     2906            }
     2907            if (upperBound(iColumn)) {
     2908              double complementarity =tVec_[iColumn]*upperSlack[iColumn];
     2909              if (complementarity>nextMu) {
     2910                double change=
     2911                  min(dualInfeasibility,
     2912                      (complementarity-nextMu)/upperSlack[iColumn]);
     2913                dualInfeasibility -= change;
     2914                printf("%d ub hicomp %g - dual inf from %g to %g\n",
     2915                       iColumn,complementarity,dualInfeasibility+change,
     2916                       dualInfeasibility);
     2917                tVec_[iColumn] -= change;
     2918                tValue = max(tVec_[iColumn],1.0e-12);
     2919              }
     2920            }
     2921          } else {
     2922            // To improve we could reduce z and/or increase t
     2923            if (lowerBound(iColumn)) {
     2924              double complementarity =zVec_[iColumn]*lowerSlack[iColumn];
     2925              if (complementarity>nextMu) {
     2926                double change=
     2927                  max(dualInfeasibility,
     2928                      (nextMu-complementarity)/lowerSlack[iColumn]);
     2929                dualInfeasibility -= change;
     2930                printf("%d lb hicomp %g - dual inf from %g to %g\n",
     2931                       iColumn,complementarity,dualInfeasibility+change,
     2932                       dualInfeasibility);
     2933                zVec_[iColumn] += change;
     2934                zValue = max(zVec_[iColumn],1.0e-12);
     2935              }
     2936            }
     2937            if (upperBound(iColumn)) {
     2938              double complementarity =tVec_[iColumn]*upperSlack[iColumn];
     2939              if (complementarity<nextMu) {
     2940                double change=
     2941                  max(dualInfeasibility,
     2942                      (complementarity-nextMu)/upperSlack[iColumn]);
     2943                dualInfeasibility -= change;
     2944                printf("%d ub locomp %g - dual inf from %g to %g\n",
     2945                       iColumn,complementarity,dualInfeasibility+change,
     2946                       dualInfeasibility);
     2947                tVec_[iColumn] -= change;
     2948                tValue = max(tVec_[iColumn],1.0e-12);
     2949              }
     2950            }
     2951          }
     2952#endif
     2953          dualFake+=newPrimal*dualInfeasibility;
     2954        }
     2955        if (lowerBoundInfeasibility>maximumBoundInfeasibility) {
     2956          maximumBoundInfeasibility=lowerBoundInfeasibility;
     2957        }
     2958        if (upperBoundInfeasibility>maximumBoundInfeasibility) {
     2959          maximumBoundInfeasibility=upperBoundInfeasibility;
     2960        }
     2961        dualInfeasibility=fabs(dualInfeasibility);
     2962        if (dualInfeasibility>maximumDualError) {
     2963          //printf("bad dual %d %g\n",iColumn,
     2964          // reducedCost-zVec[iColumn]+tVec[iColumn]+gammaTerm*newPrimal);
     2965          maximumDualError=dualInfeasibility;
     2966        }
    26082967        dualObjectiveValue+=dualObjectiveThis;
    2609         if (s>largeGap) {
    2610           s=largeGap;
    2611         }
    2612         if (t>largeGap) {
    2613           t=largeGap;
    2614         }
    2615         double divisor = s*tValue+t*zValue;
    2616         double diagonalValue=(t*s)/divisor;
     2968        gammaOffset += newPrimal*newPrimal;
     2969        if (sLower>largeGap) {
     2970          sLower=largeGap;
     2971        }
     2972        if (sUpper>largeGap) {
     2973          sUpper=largeGap;
     2974        }
     2975#if 1
     2976        double divisor = sLower*tValue+sUpper*zValue+gammaTerm*sLower*sUpper;
     2977        double diagonalValue=(sUpper*sLower)/divisor;
     2978#else
     2979        double divisor = sLower*tValue+sUpper*zValue+gammaTerm*sLower*sUpper;
     2980        double diagonalValue2=(sUpper*sLower)/divisor;
     2981        double diagonalValue;
     2982        if (!lowerBound(iColumn)) {
     2983          diagonalValue = tValue/sUpper + gammaTerm;
     2984        } else if (!upperBound(iColumn)) {
     2985          diagonalValue = zValue/sLower + gammaTerm;
     2986        } else {
     2987          diagonalValue = zValue/sLower + tValue/sUpper + gammaTerm;
     2988        }
     2989        diagonalValue = 1.0/diagonalValue;
     2990#endif
    26172991        diagonal[iColumn]=diagonalValue;
    26182992        //FUDGE
    26192993        if (diagonalValue>diagonalLimit) {
    2620           //cout<<"large diagonal "<<diagonalValue<<endl;
     2994          std::cout<<"large diagonal "<<diagonalValue<<std::endl;
    26212995          diagonal[iColumn]=diagonalLimit;
    26222996        }
    26232997        if (diagonalValue<1.0e-10) {
    2624           //cout<<"small diagonal "<<diagonalValue<<endl;
     2998          //std::cout<<"small diagonal "<<diagonalValue<<std::endl;
    26252999        }
    26263000        if (diagonalValue>largestDiagonal) {
     
    26303004          smallestDiagonal=diagonalValue;
    26313005        }
    2632         double dualInfeasibility=reducedCost-zVec[iColumn]+tVec[iColumn];
    2633         if (fabs(dualInfeasibility)>dualTolerance) {
    2634           dualFake+=newPrimal*dualInfeasibility;
    2635         }
    2636         dualInfeasibility=fabs(dualInfeasibility);
    2637         if (dualInfeasibility>maximumDualError) {
    2638           maximumDualError=dualInfeasibility;
    2639         }
    2640         deltaZ[iColumn]=0.0;
     3006        deltaX[iColumn]=0.0;
    26413007      } else {
    26423008        numberKilled++;
     
    26463012        setFlagged(iColumn);
    26473013        setFixedOrFree(iColumn);
    2648         deltaZ[iColumn]=newPrimal;
     3014        deltaX[iColumn]=newPrimal;
    26493015        offsetObjective+=newPrimal*cost[iColumn];
    26503016      }
    26513017    } else {
    2652       deltaZ[iColumn]=primal[iColumn];
     3018      deltaX[iColumn]=primal[iColumn];
    26533019      diagonal[iColumn]=0.0;
    26543020      offsetObjective+=primal[iColumn]*cost[iColumn];
     
    26663032      }
    26673033    }
    2668   }
     3034    primalObjectiveValue+=primal[iColumn]*cost[iColumn];
     3035  }
    26693036  handler_->message(CLP_BARRIER_DIAGONAL,messages_)
    26703037    <<largestDiagonal<<smallestDiagonal
     
    26813048#endif
    26823049  // update rhs region
    2683   multiplyAdd(deltaZ+numberColumns_,numberRows_,-1.0,rhsFixRegion_,1.0);
    2684   matrix_->times(1.0,deltaZ,rhsFixRegion_);
    2685   primalObjectiveValue=primalObjectiveValue+offsetObjective;
     3050  multiplyAdd(deltaX+numberColumns_,numberRows_,-1.0,rhsFixRegion_,1.0);
     3051  matrix_->times(1.0,deltaX,rhsFixRegion_);
     3052  primalObjectiveValue += 0.5*gamma2*gammaOffset+0.5*quadraticOffset;
     3053  if (quadraticOffset) {
     3054    //printf("gamma offset %g %g, quadoffset %g\n",gammaOffset,gamma2*gammaOffset,quadraticOffset);
     3055  }
     3056 
    26863057  dualObjectiveValue+=offsetObjective+dualFake;
     3058  dualObjectiveValue -= 0.5*gamma2*gammaOffset+0.5*quadraticOffset;
    26873059  if (numberIncreased||numberDecreased) {
    26883060    handler_->message(CLP_BARRIER_SLACKS,messages_)
     
    28473219  return product;
    28483220}
     3221//See exactly what would happen given current deltas
     3222void
     3223ClpPredictorCorrector::debugMove(int phase,double primalStep, double dualStep)
     3224{
     3225  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_;
     3238  int numberTotal = numberRows_+numberColumns_;
     3239  double * dualNew = ClpCopyOfArray(dual_,numberRows_);
     3240  double * errorRegionNew = new double [numberRows_];
     3241  double * rhsFixRegionNew = new double [numberRows_];
     3242  double * primalNew = ClpCopyOfArray(solution_,numberTotal);
     3243  double * djNew = new double[numberTotal];
     3244  //update pi
     3245  multiplyAdd(deltaY_,numberRows_,dualStep,dualNew,1.0);
     3246  // do reduced costs
     3247  CoinMemcpyN(dualNew,numberRows_,djNew+numberColumns_);
     3248  CoinMemcpyN(cost_,numberColumns_,djNew);
     3249  matrix_->transposeTimes(-1.0,dualNew,djNew);
     3250  double quadraticOffset=quadraticDjs(djNew,primalNew,1.0);
     3251  // update x
     3252  int iColumn;
     3253  for (iColumn=0;iColumn<numberTotal;iColumn++) {
     3254    if (!flagged(iColumn))
     3255      primalNew[iColumn] +=primalStep*deltaX[iColumn];
     3256  }
     3257  CoinZeroN(errorRegionNew,numberRows_);
     3258  CoinZeroN(rhsFixRegionNew,numberRows_);
     3259  double maximumBoundInfeasibility=0.0;
     3260  double maximumDualError=1.0e-12;
     3261  double primalObjectiveValue=0.0;
     3262  double dualObjectiveValue=0.0;
     3263  double solutionNorm=1.0e-12;
     3264  const double largeFactor=1.0e2;
     3265  double largeGap=largeFactor*solutionNorm_;
     3266  if (largeGap<largeFactor) {
     3267    largeGap=largeFactor;
     3268  }
     3269  double dualFake=0.0;
     3270  double dualTolerance =  dblParam_[ClpDualTolerance];
     3271  dualTolerance=dualTolerance/scaleFactor_;
     3272  if (dualTolerance<1.0e-12) {
     3273    dualTolerance=1.0e-12;
     3274  }
     3275  double newGap=0.0;
     3276  double offsetObjective=0.0;
     3277  double gamma2 = gamma_*gamma_; // gamma*gamma will be added to diagonal
     3278  double gammaOffset=0.0;
     3279  double maximumDjInfeasibility=0.0;
     3280  for (int iColumn=0;iColumn<numberTotal;iColumn++) {
     3281    if (!flagged(iColumn)) {
     3282      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];
     3286      double oldPrimal = solution_[iColumn];
     3287      double newPrimal=primalNew[iColumn];
     3288      double lowerBoundInfeasibility=0.0;
     3289      double upperBoundInfeasibility=0.0;
     3290      if (lowerBound(iColumn)) {
     3291        double oldSlack = lowerSlack[iColumn];
     3292        double newSlack=
     3293          lowerSlack[iColumn]+primalStep*(oldPrimal-oldSlack
     3294                                                 + thisWeight-lower[iColumn]);
     3295        if (zValue>dualTolerance) {
     3296          dualObjectiveValue+=lower[iColumn]*zVec[iColumn];
     3297        }
     3298        lowerBoundInfeasibility = fabs(newPrimal-newSlack-lower[iColumn]);
     3299        newGap += newSlack*zValue;
     3300      }
     3301      if (upperBound(iColumn)) {
     3302        double oldSlack = upperSlack[iColumn];
     3303        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;
     3311      }
     3312      if (fabs(newPrimal)>solutionNorm) {
     3313        solutionNorm=fabs(newPrimal);
     3314      }
     3315      double gammaTerm = gamma2;
     3316      if (primalR_) {
     3317        gammaTerm += primalR_[iColumn];
     3318        quadraticOffset += newPrimal*newPrimal*primalR_[iColumn];
     3319      }
     3320      double dualInfeasibility=
     3321        reducedCost-zValue+tValue+gammaTerm*newPrimal;
     3322      if (fabs(dualInfeasibility)>dualTolerance) {
     3323        dualFake+=newPrimal*dualInfeasibility;
     3324      }
     3325      if (lowerBoundInfeasibility>maximumBoundInfeasibility) {
     3326        maximumBoundInfeasibility=lowerBoundInfeasibility;
     3327      }
     3328      if (upperBoundInfeasibility>maximumBoundInfeasibility) {
     3329        maximumBoundInfeasibility=upperBoundInfeasibility;
     3330      }
     3331      dualInfeasibility=fabs(dualInfeasibility);
     3332      if (dualInfeasibility>maximumDualError) {
     3333        //printf("bad dual %d %g\n",iColumn,
     3334        // reducedCost-zVec[iColumn]+tVec[iColumn]+gammaTerm*newPrimal);
     3335        maximumDualError=dualInfeasibility;
     3336      }
     3337      gammaOffset += newPrimal*newPrimal;
     3338      djNew[iColumn]=0.0;
     3339    } 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) {
     3343          if (-djNew[iColumn]>maximumDjInfeasibility) {
     3344            maximumDjInfeasibility=-djNew[iColumn];
     3345          }
     3346        }
     3347        if (primalNew[iColumn]>upper[iColumn]-1.0e-8&&djNew[iColumn]>1.0e-8) {
     3348          if (djNew[iColumn]>maximumDjInfeasibility) {
     3349            maximumDjInfeasibility=djNew[iColumn];
     3350          }
     3351        }
     3352      }
     3353      djNew[iColumn]=primalNew[iColumn];
     3354    }
     3355    primalObjectiveValue+=primal[iColumn]*cost[iColumn];
     3356  }
     3357  // update rhs region
     3358  multiplyAdd(djNew+numberColumns_,numberRows_,-1.0,rhsFixRegionNew,1.0);
     3359  matrix_->times(1.0,djNew,rhsFixRegionNew);
     3360  primalObjectiveValue += 0.5*gamma2*gammaOffset+0.5*quadraticOffset;
     3361  dualObjectiveValue+=offsetObjective+dualFake;
     3362  dualObjectiveValue -= 0.5*gamma2*gammaOffset+0.5*quadraticOffset;
     3363  // Need to rethink (but it is only for printing)
     3364  //compute error and fixed RHS
     3365  multiplyAdd(primalNew+numberColumns_,numberRows_,-1.0,errorRegionNew,0.0);
     3366  matrix_->times(1.0,primalNew,errorRegionNew);
     3367  //finish off objective computation
     3368  double primalObjectiveNew=primalObjectiveValue*scaleFactor_;
     3369  double dualValue2=innerProduct(dualNew,numberRows_,
     3370                rhsFixRegionNew);
     3371  dualObjectiveValue-=dualValue2;
     3372  double dualObjectiveNew=dualObjectiveValue*scaleFactor_;
     3373  double maximumRHSError1=0.0;
     3374  double maximumRHSError2=0.0;
     3375  double primalOffset=0.0;
     3376  char * dropped = cholesky_->rowsDropped();
     3377  int iRow;
     3378  for (iRow=0;iRow<numberRows_;iRow++) {
     3379    double value=errorRegionNew[iRow];
     3380    if (!dropped[iRow]) {
     3381      if (fabs(value)>maximumRHSError1) {
     3382        maximumRHSError1=fabs(value);
     3383      }
     3384    } else {
     3385      if (fabs(value)>maximumRHSError2) {
     3386        maximumRHSError2=fabs(value);
     3387      }
     3388      primalOffset+=value*dualNew[iRow];
     3389    }
     3390  }
     3391  primalObjectiveNew-=primalOffset*scaleFactor_;
     3392  double maximumRHSError;
     3393  if (maximumRHSError1>maximumRHSError2) {
     3394    maximumRHSError=maximumRHSError1;
     3395  } else {
     3396    maximumRHSError=maximumRHSError1; //note change
     3397    if (maximumRHSError2>primalTolerance()) {
     3398      handler_->message(CLP_BARRIER_ABS_DROPPED,messages_)
     3399        <<maximumRHSError2
     3400        <<CoinMessageEol;
     3401    }
     3402  }
     3403  printf("PH %d %g, %g new comp %g, b %g, p %g, d %g\n",phase,
     3404         primalStep,dualStep,newGap,maximumBoundInfeasibility,
     3405         maximumRHSError,maximumDualError);
     3406  if (handler_->logLevel()>1)
     3407    printf("       objs %g %g\n",
     3408           primalObjectiveNew,dualObjectiveNew);
     3409  if (maximumDjInfeasibility) {
     3410    printf(" max dj error on fixed %g\n",
     3411           maximumDjInfeasibility);
     3412  }
     3413  delete [] dualNew;
     3414  delete [] errorRegionNew;
     3415  delete [] rhsFixRegionNew;
     3416  delete [] primalNew;
     3417  delete [] djNew;
     3418}
  • trunk/ClpPrimalColumnSteepest.cpp

    r341 r389  
    27772777  if (mode==1&&weights_) {
    27782778    if (pivotSequence_>=0) {
     2779      // clean array first
     2780      alternateWeights_->clear();
    27792781      // save pivot order
    27802782      memcpy(alternateWeights_->getIndices(),pivotVariable,
     
    28452847        pivotSequence_= -1;
    28462848        model_->setSequenceOut(-1);
     2849        alternateWeights_->clear();
    28472850      }
    28482851    }
     
    28552858  }
    28562859  if (mode>=2&&mode!=5) {
    2857     if (mode!=3&&pivotSequence_>=0) {
    2858       // restore pivot row
    2859       int iRow;
    2860       // permute alternateWeights
    2861       double * temp = new double[numberRows+numberColumns];
    2862       double * work = alternateWeights_->denseVector();
    2863       int * oldPivotOrder = alternateWeights_->getIndices();
    2864       for (iRow=0;iRow<numberRows;iRow++) {
    2865         int iPivot=oldPivotOrder[iRow];
    2866         temp[iPivot]=work[iRow];
    2867       }
    2868       int number=0;
    2869       int found=-1;
    2870       int * which = oldPivotOrder;
    2871       // find pivot row and re-create alternateWeights
    2872       for (iRow=0;iRow<numberRows;iRow++) {
    2873         int iPivot=pivotVariable[iRow];
    2874         if (iPivot==pivotSequence_)
    2875           found=iRow;
    2876         work[iRow]=temp[iPivot];
    2877         if (work[iRow])
    2878           which[number++]=iRow;
    2879       }
    2880       alternateWeights_->setNumElements(number);
     2860    if (mode!=3) {
     2861      if (pivotSequence_>=0) {
     2862        // restore pivot row
     2863        int iRow;
     2864        // permute alternateWeights
     2865        double * temp = new double[numberRows+numberColumns];
     2866        double * work = alternateWeights_->denseVector();
     2867        int * oldPivotOrder = alternateWeights_->getIndices();
     2868        for (iRow=0;iRow<numberRows;iRow++) {
     2869          int iPivot=oldPivotOrder[iRow];
     2870          temp[iPivot]=work[iRow];
     2871        }
     2872        int number=0;
     2873        int found=-1;
     2874        int * which = oldPivotOrder;
     2875        // find pivot row and re-create alternateWeights
     2876        for (iRow=0;iRow<numberRows;iRow++) {
     2877          int iPivot=pivotVariable[iRow];
     2878          if (iPivot==pivotSequence_)
     2879            found=iRow;
     2880          work[iRow]=temp[iPivot];
     2881          if (work[iRow])
     2882            which[number++]=iRow;
     2883        }
     2884        alternateWeights_->setNumElements(number);
    28812885#ifdef CLP_DEBUG
    2882       // Can happen but I should clean up
    2883       assert(found>=0);
     2886        // Can happen but I should clean up
     2887        assert(found>=0);
    28842888#endif
    2885       pivotSequence_ = found;
    2886       delete [] temp;
     2889        pivotSequence_ = found;
     2890        delete [] temp;
     2891      } else {
     2892        // Just clean up
     2893        alternateWeights_->clear();
     2894      }
    28872895    }
    28882896    // Save size of factorization
  • trunk/ClpQuadraticObjective.cpp

    r277 r389  
    33
    44#include "CoinPragma.hpp"
     5#include "CoinHelperFunctions.hpp"
    56#include "ClpModel.hpp"
    67#include "ClpQuadraticObjective.hpp"
    7 
    88//#############################################################################
    99// Constructors / Destructor / Assignment
     
    5959// Copy constructor
    6060//-------------------------------------------------------------------
    61 ClpQuadraticObjective::ClpQuadraticObjective (const ClpQuadraticObjective & rhs)
     61ClpQuadraticObjective::ClpQuadraticObjective (const ClpQuadraticObjective & rhs,
     62                                              int type)
    6263: ClpObjective(rhs)
    6364
     
    7677    gradient_=NULL;
    7778  }
    78   if (rhs.quadraticObjective_)
    79     quadraticObjective_ = new CoinPackedMatrix(*rhs.quadraticObjective_);
    80   else
     79  if (rhs.quadraticObjective_) {
     80    // see what type of matrix wanted
     81    if (type==0) {
     82      // just copy
     83      quadraticObjective_ = new CoinPackedMatrix(*rhs.quadraticObjective_);
     84    } else if (type==1) {
     85      // expand to full symmetric
     86      const int * columnQuadratic1 = rhs.quadraticObjective_->getIndices();
     87      const CoinBigIndex * columnQuadraticStart1 = rhs.quadraticObjective_->getVectorStarts();
     88      const int * columnQuadraticLength1 = rhs.quadraticObjective_->getVectorLengths();
     89      const double * quadraticElement1 = rhs.quadraticObjective_->getElements();
     90      CoinBigIndex * columnQuadraticStart2 = new CoinBigIndex [numberExtendedColumns_+1];
     91      int * columnQuadraticLength2 = new int [numberExtendedColumns_];
     92      int iColumn;
     93      int numberColumns = rhs.quadraticObjective_->getNumCols();
     94      int numberBelow=0;
     95      int numberAbove=0;
     96      int numberDiagonal=0;
     97      CoinZeroN(columnQuadraticLength2,numberExtendedColumns_);
     98      for (iColumn=0;iColumn<numberColumns;iColumn++) {
     99        for (CoinBigIndex j=columnQuadraticStart1[iColumn];
     100             j<columnQuadraticStart1[iColumn]+columnQuadraticLength1[iColumn];j++) {
     101          int jColumn = columnQuadratic1[j];
     102          if (jColumn>iColumn) {
     103            numberBelow++;
     104            columnQuadraticLength2[jColumn]++;
     105            columnQuadraticLength2[iColumn]++;
     106          } else if (jColumn==iColumn) {
     107            numberDiagonal++;
     108            columnQuadraticLength2[iColumn]++;
     109          } else {
     110            numberAbove++;
     111          }
     112        }
     113      }
     114      if (numberAbove>0) {
     115        if (numberAbove==numberBelow) {
     116          // already done
     117          quadraticObjective_ = new CoinPackedMatrix(*rhs.quadraticObjective_);
     118          delete [] columnQuadraticStart2;
     119          delete [] columnQuadraticLength2;
     120        } else {
     121          printf("number above = %d, number below = %d, error\n",
     122                 numberAbove,numberBelow);
     123        }
     124      } else {
     125        int numberElements=numberDiagonal+2*numberBelow;
     126        int * columnQuadratic2 = new int [numberElements];
     127        double * quadraticElement2 = new double [numberElements];
     128        columnQuadraticStart2[0]=0;
     129        numberElements=0;
     130        for (iColumn=0;iColumn<numberColumns;iColumn++) {
     131          int n=columnQuadraticLength2[iColumn];
     132          columnQuadraticLength2[iColumn]=0;
     133          numberElements += n;
     134          columnQuadraticStart2[iColumn+1]=numberElements;
     135        }
     136        for (iColumn=0;iColumn<numberColumns;iColumn++) {
     137          for (CoinBigIndex j=columnQuadraticStart1[iColumn];
     138               j<columnQuadraticStart1[iColumn]+columnQuadraticLength1[iColumn];j++) {
     139            int jColumn = columnQuadratic1[j];
     140            if (jColumn>iColumn) {
     141              // put in two places
     142              CoinBigIndex put=columnQuadraticLength2[jColumn]+columnQuadraticStart2[jColumn];
     143              columnQuadraticLength2[jColumn]++;
     144              quadraticElement2[put]=quadraticElement1[j];
     145              columnQuadratic2[put]=iColumn;
     146              put=columnQuadraticLength2[iColumn]+columnQuadraticStart2[iColumn];
     147              columnQuadraticLength2[iColumn]++;
     148              quadraticElement2[put]=quadraticElement1[j];
     149              columnQuadratic2[put]=jColumn;
     150            } else if (jColumn==iColumn) {
     151              CoinBigIndex put=columnQuadraticLength2[iColumn]+columnQuadraticStart2[iColumn];
     152              columnQuadraticLength2[iColumn]++;
     153              quadraticElement2[put]=quadraticElement1[j];
     154              columnQuadratic2[put]=iColumn;
     155            } else {
     156              abort();
     157            }
     158          }
     159        }
     160        // Now create
     161        quadraticObjective_ =
     162          new CoinPackedMatrix (true,
     163                                rhs.numberExtendedColumns_,
     164                                rhs.numberExtendedColumns_,
     165                                numberElements,
     166                                quadraticElement2,
     167                                columnQuadratic2,
     168                                columnQuadraticStart2,
     169                                columnQuadraticLength2,0.0,0.0);
     170        delete [] columnQuadraticStart2;
     171        delete [] columnQuadraticLength2;
     172        delete [] columnQuadratic2;
     173        delete [] quadraticElement2;
     174      }
     175    } else {
     176      abort(); // code when needed
     177    }
     178           
     179  } else {
    81180    quadraticObjective_=NULL;
     181  }
    82182}
    83183/* Subset constructor.  Duplicates are allowed
  • trunk/ClpSimplex.cpp

    r383 r389  
    4444  remainingDualInfeasibility_(0.0),
    4545  largeValue_(1.0e15),
    46   objectiveScale_(1.0),
    47   rhsScale_(1.0),
    4846  largestPrimalError_(0.0),
    4947  largestDualError_(0.0),
     
    149147  remainingDualInfeasibility_(0.0),
    150148  largeValue_(1.0e15),
    151   objectiveScale_(1.0),
    152   rhsScale_(1.0),
    153149  largestPrimalError_(0.0),
    154150  largestDualError_(0.0),
     
    13641360  remainingDualInfeasibility_(0.0),
    13651361  largeValue_(1.0e15),
    1366   objectiveScale_(1.0),
    1367   rhsScale_(1.0),
    13681362  largestPrimalError_(0.0),
    13691363  largestDualError_(0.0),
     
    14631457  remainingDualInfeasibility_(0.0),
    14641458  largeValue_(1.0e15),
    1465   objectiveScale_(1.0),
    1466   rhsScale_(1.0),
    14671459  largestPrimalError_(0.0),
    14681460  largestDualError_(0.0),
     
    16241616  remainingDualInfeasibility_ = rhs.remainingDualInfeasibility_;
    16251617  largeValue_ = rhs.largeValue_;
    1626   objectiveScale_ = rhs.objectiveScale_;
    1627   rhsScale_ = rhs.rhsScale_;
    16281618  largestPrimalError_ = rhs.largestPrimalError_;
    16291619  largestDualError_ = rhs.largestDualError_;
     
    34733463  return 0;
    34743464}
     3465/* Primal ranging.
     3466   This computes increase/decrease in value for each given variable and corresponding
     3467   sequence numbers which would change basis.  Sequence numbers are 0..numberColumns
     3468   and numberColumns.. for artificials/slacks.
     3469   For basic variables the sequence number will be that of the basic variables.
     3470   
     3471   Up to user to providen correct length arrays.
     3472   
     3473   Returns non-zero if infeasible unbounded etc
     3474*/
     3475int ClpSimplex::primalRanging(int numberCheck,const int * which,
     3476                  double * valueIncrease, int * sequenceIncrease,
     3477                  double * valueDecrease, int * sequenceDecrease)
     3478{
     3479  int savePerturbation = perturbation_;
     3480  perturbation_=100;
     3481  int returnCode = ((ClpSimplexPrimal *) this)->primal(0,1);
     3482  if (problemStatus_==10) {
     3483    //printf("Cleaning up with dual\n");
     3484    bool denseFactorization = initialDenseFactorization();
     3485    // It will be safe to allow dense
     3486    setInitialDenseFactorization(true);
     3487    // check which algorithms allowed
     3488    int dummy;
     3489    if ((matrix_->generalExpanded(this,4,dummy)&2)!=0) {
     3490      // upperOut_ has largest away from bound
     3491      double saveBound=dualBound_;
     3492      if (upperOut_>0.0)
     3493        dualBound_=2.0*upperOut_;
     3494      returnCode = ((ClpSimplexDual *) this)->dual(0,1);
     3495      dualBound_=saveBound;
     3496    } else {
     3497      returnCode = ((ClpSimplexPrimal *) this)->primal(0,1);
     3498    }
     3499    setInitialDenseFactorization(denseFactorization);
     3500    if (problemStatus_==10)
     3501      problemStatus_=0;
     3502  }
     3503  perturbation_=savePerturbation;
     3504  if (problemStatus_) {
     3505    finish(); // get rid of arrays
     3506    return 1; // odd status
     3507  }
     3508  ((ClpSimplexOther *) this)->primalRanging(numberCheck,which,
     3509                                          valueIncrease,sequenceIncrease,
     3510                                          valueDecrease,sequenceDecrease);
     3511  finish(); // get rid of arrays
     3512  return 0;
     3513}
    34753514#include "ClpSimplexPrimalQuadratic.hpp"
    34763515/* Solves quadratic problem using SLP - may be used as crash
     
    35293568  ClpSimplex * saveModel2=NULL;
    35303569  int numberFixed = barrier.numberFixed();
    3531   if (numberFixed*20>barrier.numberRows()&&numberFixed>5000&&crossover) {
     3570  if (numberFixed*20>barrier.numberRows()&&numberFixed>5000&&crossover&&0) {
    35323571    // may as well do presolve
    35333572    int numberRows = barrier.numberRows();
     
    36093648      int numberColumns = model2->numberColumns();
    36103649      // just primal values pass
     3650      double saveScale = model2->objectiveScale();
     3651      model2->setObjectiveScale(1.0e-3);
    36113652      model2->primal(2);
     3653      model2->setObjectiveScale(saveScale);
    36123654      // save primal solution and copy back dual
    36133655      CoinMemcpyN(model2->primalRowSolution(),
     
    36793721                  numberColumns,model2->primalColumnSolution());
    36803722    }
     3723//     double saveScale = model2->objectiveScale();
     3724//     model2->setObjectiveScale(1.0e-3);
     3725//     model2->primal(2);
     3726//    model2->setObjectiveScale(saveScale);
    36813727    model2->primal(1);
    36823728  } else if (barrierStatus==4&&crossover) {
     
    37553801  specialOptions_ = otherModel.specialOptions_;
    37563802  automaticScale_ = otherModel.automaticScale_;
    3757   objectiveScale_ = otherModel.objectiveScale_;
    3758   rhsScale_ = otherModel.rhsScale_;
    37593803}
    37603804typedef struct {
     
    54415485    objective_[i] = COIN_DBL_MAX;
    54425486    infeasibility_[i] = -1.0; // set to an impossible value
     5487    realInfeasibility_[i] = COIN_DBL_MAX;
    54435488    numberInfeasibilities_[i]=-1;
    54445489    iterationNumber_[i]=-1;
     
    54695514    objective_[i] = rhs.objective_[i];
    54705515    infeasibility_[i] = rhs.infeasibility_[i];
     5516    realInfeasibility_[i] = rhs.realInfeasibility_[i];
    54715517    numberInfeasibilities_[i]=rhs.numberInfeasibilities_[i];
    54725518    iterationNumber_[i]=rhs.iterationNumber_[i];
     
    54945540      objective_[i] = -COIN_DBL_MAX;
    54955541    infeasibility_[i] = -1.0; // set to an impossible value
     5542    realInfeasibility_[i] = COIN_DBL_MAX;
    54965543    numberInfeasibilities_[i]=-1;
    54975544    iterationNumber_[i]=-1;
     
    55165563      objective_[i] = rhs.objective_[i];
    55175564      infeasibility_[i] = rhs.infeasibility_[i];
     5565      realInfeasibility_[i] = rhs.realInfeasibility_[i];
    55185566      numberInfeasibilities_[i]=rhs.numberInfeasibilities_[i];
    55195567      iterationNumber_[i]=rhs.iterationNumber_[i];
     
    55495597  double objective = model_->objectiveValue();
    55505598  double infeasibility;
     5599  double realInfeasibility=0.0;
    55515600  int numberInfeasibilities;
    55525601  int iterationNumber = model_->numberIterations();
     
    55585607    //primal
    55595608    infeasibility = model_->sumDualInfeasibilities();
     5609    realInfeasibility = model_->nonLinearCost()->sumInfeasibilities();
    55605610    numberInfeasibilities = model_->numberDualInfeasibilities();
    55615611  }
     
    55885638      objective_[i-1] = objective_[i];
    55895639      infeasibility_[i-1] = infeasibility_[i];
     5640      realInfeasibility_[i-1] = realInfeasibility_[i];
    55905641      numberInfeasibilities_[i-1]=numberInfeasibilities_[i];
    55915642      iterationNumber_[i-1]=iterationNumber_[i];
     
    55945645  objective_[CLP_PROGRESS-1] = objective;
    55955646  infeasibility_[CLP_PROGRESS-1] = infeasibility;
     5647  realInfeasibility_[CLP_PROGRESS-1] = realInfeasibility;
    55965648  numberInfeasibilities_[CLP_PROGRESS-1]=numberInfeasibilities;
    55975649  iterationNumber_[CLP_PROGRESS-1]=iterationNumber;
     
    56905742  return objective_[CLP_PROGRESS-1-back];
    56915743}
     5744// Returns previous infeasibility (if -1) - current if (0)
     5745double
     5746ClpSimplexProgress::lastInfeasibility(int back) const
     5747{
     5748  return realInfeasibility_[CLP_PROGRESS-1-back];
     5749}
     5750// Sets real primal infeasibility
     5751void
     5752ClpSimplexProgress::setInfeasibility(double value)
     5753{
     5754  for (int i=1;i<CLP_PROGRESS;i++)
     5755    realInfeasibility_[i-1] = realInfeasibility_[i];
     5756  realInfeasibility_[CLP_PROGRESS-1]=value;
     5757}
    56925758// Modify objective e.g. if dual infeasible in dual
    56935759void
  • trunk/ClpSimplexOther.cpp

    r386 r389  
    6060                                rowArray_[0],rowArray_[3],columnArray_[0]);
    6161        // do ratio test up and down
    62         checkRatios(rowArray_[0],columnArray_[0],costIncrease,sequenceIncrease,
     62        checkDualRatios(rowArray_[0],columnArray_[0],costIncrease,sequenceIncrease,
    6363                    costDecrease,sequenceDecrease);
    6464      }
     
    124124*/
    125125void
    126 ClpSimplexOther::checkRatios(CoinIndexedVector * rowArray,
     126ClpSimplexOther::checkDualRatios(CoinIndexedVector * rowArray,
    127127                             CoinIndexedVector * columnArray,
    128128                             double & costIncrease, int & sequenceIncrease,
     
    221221  }
    222222}
     223/** Primal ranging.
     224    This computes increase/decrease in value for each given variable and corresponding
     225    sequence numbers which would change basis.  Sequence numbers are 0..numberColumns
     226    and numberColumns.. for artificials/slacks.
     227    For basic variables the sequence number will be that of the basic variables.
     228   
     229    Up to user to provide correct length arrays.
     230   
     231    When here - guaranteed optimal
     232*/
     233void
     234ClpSimplexOther::primalRanging(int numberCheck,const int * which,
     235                  double * valueIncreased, int * sequenceIncreased,
     236                  double * valueDecreased, int * sequenceDecreased)
     237{
     238  rowArray_[0]->clear();
     239  rowArray_[1]->clear();
     240  lowerIn_=-COIN_DBL_MAX;
     241  upperIn_=COIN_DBL_MAX;
     242  valueIn_ = 0.0;
     243  for ( int i=0;i<numberCheck;i++) {
     244    int iSequence = which[i];
     245    double valueIncrease=COIN_DBL_MAX;
     246    double valueDecrease=COIN_DBL_MAX;
     247    int sequenceIncrease=-1;
     248    int sequenceDecrease=-1;
     249   
     250    switch(getStatus(iSequence)) {
     251     
     252    case basic:
     253    case isFree:
     254    case superBasic:
     255      // Easy
     256      valueIncrease=max(0.0,upper_[iSequence]-solution_[iSequence]);
     257      valueDecrease=max(0.0,solution_[iSequence]-lower_[iSequence]);
     258      sequenceIncrease=iSequence;
     259      sequenceDecrease=iSequence;
     260      break;
     261    case isFixed:
     262    case atUpperBound:
     263    case atLowerBound:
     264      {
     265        // Non trivial
     266        // Other bound is ignored
     267        unpackPacked(rowArray_[1],iSequence);
     268        factorization_->updateColumn(rowArray_[2],rowArray_[1]);
     269        // Get extra rows
     270        matrix_->extendUpdated(this,rowArray_[1],0);
     271        // do ratio test
     272        checkPrimalRatios(rowArray_[1],1);
     273        if (pivotRow_>=0) {
     274          valueIncrease = theta_;
     275          sequenceIncrease=pivotVariable_[pivotRow_];
     276        }
     277        directionIn_=-1; // down
     278        checkPrimalRatios(rowArray_[1],-1);
     279        if (pivotRow_>=0) {
     280          valueDecrease = theta_;
     281          sequenceDecrease=pivotVariable_[pivotRow_];
     282        }
     283        rowArray_[1]->clear();
     284      }
     285      break;
     286    }
     287    double scaleFactor;
     288    if (rowScale_) {
     289      if (iSequence<numberColumns_)
     290        scaleFactor = columnScale_[iSequence]/rhsScale_;
     291      else
     292        scaleFactor = 1.0/(rowScale_[iSequence-numberColumns_]*rhsScale_);
     293    } else {
     294      scaleFactor = 1.0/rhsScale_;
     295    }
     296    if (valueIncrease<1.0e30)
     297      valueIncrease *= scaleFactor;
     298    else
     299      valueIncrease = COIN_DBL_MAX;
     300    if (valueDecrease<1.0e30)
     301      valueDecrease *= scaleFactor;
     302    else
     303      valueDecrease = COIN_DBL_MAX;
     304    valueIncreased[i] = valueIncrease;
     305    sequenceIncreased[i] = sequenceIncrease;
     306    valueDecreased[i] = valueDecrease;
     307    sequenceDecreased[i] = sequenceDecrease;
     308  }
     309}
     310/*
     311   Row array has pivot column
     312   This is used in primal ranging
     313*/
     314void
     315ClpSimplexOther::checkPrimalRatios(CoinIndexedVector * rowArray,
     316                         int direction)
     317{
     318  // sequence stays as row number until end
     319  pivotRow_=-1;
     320  double acceptablePivot=1.0e-7;
     321  double * work=rowArray->denseVector();
     322  int number=rowArray->getNumElements();
     323  int * which=rowArray->getIndices();
     324
     325  // we need to swap sign if going down
     326  double way = direction;
     327  theta_ = 1.0e30;
     328  for (int iIndex=0;iIndex<number;iIndex++) {
     329   
     330    int iRow = which[iIndex];
     331    double alpha = work[iIndex]*way;
     332    int iPivot=pivotVariable_[iRow];
     333    double oldValue = solution_[iPivot];
     334    if (fabs(alpha)>acceptablePivot) {
     335      if (alpha>0.0) {
     336        // basic variable going towards lower bound
     337        double bound = lower_[iPivot];
     338        oldValue -= bound;
     339        if (oldValue-theta_*alpha<0.0) {
     340          pivotRow_ = iRow;
     341          theta_ = max(0.0,oldValue/alpha);
     342        }
     343      } else {
     344        // basic variable going towards upper bound
     345        double bound = upper_[iPivot];
     346        oldValue = oldValue-bound;
     347        if (oldValue-theta_*alpha>0.0) {
     348          pivotRow_ = iRow;
     349          theta_ = max(0.0,oldValue/alpha);
     350        }
     351      }
     352    }
     353  }
     354}
  • trunk/ClpSimplexPrimal.cpp

    r372 r389  
    640640            return;
    641641          }
    642 #if 1
    643           // switch off dense
    644           int saveDense = factorization_->denseThreshold();
    645           factorization_->setDenseThreshold(0);
    646           // make sure will do safe factorization
    647           pivotVariable_[0]=-1;
    648           internalFactorize(2);
    649           factorization_->setDenseThreshold(saveDense);
    650           // restore extra stuff
    651           matrix_->generalExpanded(this,6,dummy);
    652 #else
    653          
    654           // no - restore previous basis
    655           assert (type==1);
    656           memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
    657           memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
    658                  numberRows_*sizeof(double));
    659           memcpy(columnActivityWork_,savedSolution_ ,
    660                  numberColumns_*sizeof(double));
    661           // restore extra stuff
    662           matrix_->generalExpanded(this,6,dummy);
    663           matrix_->generalExpanded(this,5,dummy);
    664           forceFactorization_=1; // a bit drastic but ..
    665           type = 2;
    666           assert (internalFactorize(1)==0);
    667 #endif
     642          if (type!=1||largestPrimalError_>1.0e3
     643              ||largestDualError_>1.0e3) {
     644            // switch off dense
     645            int saveDense = factorization_->denseThreshold();
     646            factorization_->setDenseThreshold(0);
     647            // make sure will do safe factorization
     648            pivotVariable_[0]=-1;
     649            internalFactorize(2);
     650            factorization_->setDenseThreshold(saveDense);
     651            // restore extra stuff
     652            matrix_->generalExpanded(this,6,dummy);
     653          } else {
     654            // no - restore previous basis
     655            memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
     656            memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
     657                   numberRows_*sizeof(double));
     658            memcpy(columnActivityWork_,savedSolution_ ,
     659                   numberColumns_*sizeof(double));
     660            // restore extra stuff
     661            matrix_->generalExpanded(this,6,dummy);
     662            matrix_->generalExpanded(this,5,dummy);
     663            forceFactorization_=1; // a bit drastic but ..
     664            type = 2;
     665            assert (internalFactorize(1)==0);
     666          }
    668667          changeMade_++; // say change made
    669668        }
     
    747746    nonLinearCost_->checkInfeasibilities(primalTolerance_);
    748747    gutsOfSolution(NULL,NULL,ifValuesPass!=0);
     748    nonLinearCost_->checkInfeasibilities(primalTolerance_);
     749  }
     750  double trueInfeasibility =nonLinearCost_->sumInfeasibilities();
     751  if (trueInfeasibility>1.0) {
     752    // If infeasibility going up may change weights
     753    double testValue = trueInfeasibility-1.0e-4*(10.0+trueInfeasibility);
     754    if(progress->lastInfeasibility()<testValue) {
     755      if (infeasibilityCost_<1.0e14) {
     756        infeasibilityCost_ *= 1.5;
     757        printf("increasing weight to %g\n",infeasibilityCost_);
     758        gutsOfSolution(NULL,NULL,ifValuesPass!=0);
     759      }
     760    }
    749761  }
    750762  // we may wish to say it is optimal even if infeasible
  • trunk/ClpSolve.cpp

    r385 r389  
    1313#include "CoinSort.hpp"
    1414#include "ClpFactorization.hpp"
     15#include "ClpQuadraticObjective.hpp"
    1516#include "ClpSimplex.hpp"
    1617#include "ClpInterior.hpp"
     
    100101  int doCrash=0;
    101102  int doSprint=0;
    102   if (method!=ClpSolve::useDual) {
     103  if (method!=ClpSolve::useDual&&method!=ClpSolve::useBarrier) {
    103104    switch (options.getSpecialOption(1)) {
    104105    case 0:
     
    907908    if (interrupt)
    908909      currentModel2 = &barrier;
     910    int barrierOptions = options.getSpecialOption(4);
     911    bool aggressiveGamma=false;
     912    bool scale=false;
     913    if (barrierOptions&8) {
     914      barrierOptions &= ~8;
     915      aggressiveGamma=true;
     916    }
     917    if (barrierOptions&4) {
     918      barrierOptions &= ~4;
     919      scale=true;
     920    }
    909921#ifdef REAL_BARRIER
    910     // uncomment this if you have Anshul Gupta's wsmp package
    911     ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp(max(100,model2->numberRows()/10));
    912     //ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp();
    913     //ClpCholeskyWssmpKKT * cholesky = new ClpCholeskyWssmpKKT(max(100,model2->numberRows()/10));
     922    // If quadratic force KKT
     923#ifndef NO_RTTI
     924    ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(barrier.objectiveAsObject()));
     925#else
     926    ClpQuadraticObjective * quadraticObj = NULL;
     927    if (objective_->type()==2)
     928      quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
     929#endif
     930    if (quadraticObj) {
     931      barrierOptions = 2;
     932    }
     933    switch (barrierOptions) {
     934    case 0:
     935      {
     936        ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp();
     937        barrier.setCholesky(cholesky);
     938      }
     939      break;
     940    case 1:
     941      {
     942        ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp(max(100,model2->numberRows()/10));
     943        barrier.setCholesky(cholesky);
     944      }
     945      break;
     946    case 2:
     947      {
     948        ClpCholeskyWssmpKKT * cholesky = new ClpCholeskyWssmpKKT(max(100,model2->numberRows()/10));
     949        barrier.setCholesky(cholesky);
     950      }
     951      break;
     952    case 3:
     953      break;
     954    default:
     955      abort();
     956    }
    914957    // uncomment this if you have Sivan Toledo's Taucs package
    915958    //ClpCholeskyTaucs * cholesky = new ClpCholeskyTaucs();
    916     barrier.setCholesky(cholesky);
     959    //barrier.setCholesky(cholesky);
    917960#endif
    918961    int numberRows = model2->numberRows();
     
    924967    }
    925968#ifndef SAVEIT
     969    //barrier.setGamma(1.0e-8);
     970    //barrier.setDelta(1.0e-8);
     971    if (aggressiveGamma) {
     972      barrier.setGamma(1.0e-3);
     973      barrier.setDelta(1.0e-3);
     974    }
     975    if (scale)
     976      barrier.scaling(1);
     977    else
     978      barrier.scaling(0);
    926979    barrier.primalDual();
    927     //printf("********** Stopping as this is debug run\n");
    928     //exit(99);
    929980#elif SAVEIT==1
    930981    barrier.primalDual();
     
    9561007    ClpSimplex * saveModel2=NULL;
    9571008    int numberFixed = barrier.numberFixed();
    958     if (numberFixed*20>barrier.numberRows()&&numberFixed>5000) {
     1009    if (numberFixed*20>barrier.numberRows()&&numberFixed>5000&&0) {
    9591010      // may as well do presolve
    9601011      int numberRows = barrier.numberRows();
     
    9691020      barrier.fixFixed();
    9701021      saveModel2=model2;
     1022    } else if (numberFixed) {
     1023      // Set fixed to bounds (may have restored earlier solution)
     1024      barrier.fixFixed(false);
    9711025    }
    9721026#ifdef BORROW   
     
    10131067      }
    10141068    }
    1015     if (maxIts&&barrierStatus<4) {
    1016       printf("***** crossover - needs more thought on difficult models\n");
     1069    if (maxIts&&barrierStatus<4&&!quadraticObj) {
     1070      //printf("***** crossover - needs more thought on difficult models\n");
    10171071#if SAVEIT==1
    10181072      model2->ClpSimplex::saveModel("xx.save");
     
    10701124        int numberColumns = model2->numberColumns();
    10711125        // just primal values pass
     1126        double saveScale = model2->objectiveScale();
     1127        model2->setObjectiveScale(1.0e-3);
    10721128        model2->primal(2);
     1129        model2->setObjectiveScale(saveScale);
    10731130        // save primal solution and copy back dual
    10741131        CoinMemcpyN(model2->primalRowSolution(),
     
    11401197                    numberColumns,model2->primalColumnSolution());
    11411198      }
     1199      double saveScale = model2->objectiveScale();
     1200      model2->setObjectiveScale(1.0e-3);
     1201      model2->primal(2);
     1202      model2->setObjectiveScale(saveScale);
    11421203      model2->primal(1);
    11431204#else
  • trunk/Test/ClpMain.cpp

    r383 r389  
    6363  DIRECTION=201,DUALPIVOT,SCALING,ERRORSALLOWED,KEEPNAMES,SPARSEFACTOR,
    6464  PRIMALPIVOT,PRESOLVE,CRASH,BIASLU,PERTURBATION,MESSAGES,AUTOSCALE,
     65  CHOLESKY,BARRIERSCALE,GAMMA,
    6566 
    6667  DIRECTORY=301,IMPORT,EXPORT,RESTORE,SAVE,DUALSIMPLEX,PRIMALSIMPLEX,
     
    935936    parameters[numberParameters-1].append("LL");
    936937    parameters[numberParameters++]=
     938      ClpItem("bscale","Whether to scale in barrier",
     939              "off",BARRIERSCALE,false);
     940    parameters[numberParameters-1].append("on");
     941    parameters[numberParameters++]=
     942      ClpItem("chol!esky","Which cholesky algorithm",
     943              "wssmp",CHOLESKY,false);
     944    parameters[numberParameters-1].append("fudge!Long");
     945    parameters[numberParameters-1].append("KKT");
     946    parameters[numberParameters-1].append("dense");
     947    parameters[numberParameters++]=
    937948      ClpItem("crash","Whether to create basis for problem",
    938949              "off",CRASH);
     
    10521063              1.0,1.0e15,FAKEBOUND,false);
    10531064    parameters[numberParameters++]=
     1065      ClpItem("gamma","Whether to regularize barrier",
     1066              "off",GAMMA,false);
     1067    parameters[numberParameters-1].append("on");
     1068    parameters[numberParameters++]=
    10541069      ClpItem("help","Print out version, non-standard options and some help",
    10551070              HELP);
     
    14161431    int numberGoodCommands=0;
    14171432    bool * goodModels = new bool[1];
    1418    
     1433
     1434    // Hidden stuff for barrier
     1435    int choleskyType = 0;
     1436    int gamma=0;
     1437    int scaleBarrier=0;
    14191438   
    14201439    int iModel=0;
     
    16921711              models[iModel].messageHandler()->setPrefix(action!=0);
    16931712              break;
     1713            case CHOLESKY:
     1714              choleskyType = action;
     1715              break;
     1716            case GAMMA:
     1717              gamma=action;
     1718              break;
     1719            case BARRIERSCALE:
     1720              scaleBarrier=action;
     1721              break;
    16941722            default:
    16951723              abort();
     
    17511779                  }
    17521780                }
     1781              } else if (method==ClpSolve::useBarrier) {
     1782                int barrierOptions = choleskyType;
     1783                if (scaleBarrier)
     1784                  barrierOptions |= 4;
     1785                if (gamma)
     1786                  barrierOptions |= 8;
     1787                solveOptions.setSpecialOption(4,barrierOptions);
    17531788              }
    17541789              model2->initialSolve(solveOptions);
  • trunk/Test/unitTest.cpp

    r383 r389  
    1818#include "ClpFactorization.hpp"
    1919#include "ClpSimplex.hpp"
    20 //#include "ClpInterior.hpp"
     20#include "ClpInterior.hpp"
    2121#include "ClpLinearObjective.hpp"
    2222#include "ClpDualRowSteepest.hpp"
     
    567567                         m.getRowLower(),m.getRowUpper());
    568568    model.primal();
    569     int which[8] = {0,1,2,3,4,5,6,7};
    570     double costIncrease[8];
    571     int sequenceIncrease[8];
    572     double costDecrease[8];
    573     int sequenceDecrease[8];
     569    int which[13] = {0,1,2,3,4,5,6,7,8,9,10,11,12};
     570    double costIncrease[13];
     571    int sequenceIncrease[13];
     572    double costDecrease[13];
     573    int sequenceDecrease[13];
    574574    // ranging
    575     model.dualRanging(8,which,costIncrease,sequenceIncrease,
     575    model.dualRanging(13,which,costIncrease,sequenceIncrease,
    576576                      costDecrease,sequenceDecrease);
    577577    int i;
    578     for ( i=0;i<8;i++)
     578    for ( i=0;i<13;i++)
    579579      printf("%d increase %g %d, decrease %g %d\n",
    580580             i,costIncrease[i],sequenceIncrease[i],
     
    590590        obj[j] *= -1.0;
    591591    }
    592     double costIncrease2[8];
    593     int sequenceIncrease2[8];
    594     double costDecrease2[8];
    595     int sequenceDecrease2[8];
     592    double costIncrease2[13];
     593    int sequenceIncrease2[13];
     594    double costDecrease2[13];
     595    int sequenceDecrease2[13];
    596596    // ranging
    597     model.dualRanging(8,which,costIncrease2,sequenceIncrease2,
     597    model.dualRanging(13,which,costIncrease2,sequenceIncrease2,
    598598                      costDecrease2,sequenceDecrease2);
    599     for (i=0;i<8;i++) {
     599    for (i=0;i<13;i++) {
    600600      assert (fabs(costIncrease[i]-costDecrease2[i])<1.0e-6);
    601601      assert (fabs(costDecrease[i]-costIncrease2[i])<1.0e-6);
     
    603603      assert (sequenceDecrease[i]==sequenceIncrease2[i]);
    604604    }
     605  }
     606  // Test primal ranging
     607  {   
     608    CoinMpsIO m;
     609    std::string fn = mpsDir+"exmip1";
     610    m.readMps(fn.c_str(),"mps");
     611    ClpSimplex model;
     612    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
     613                         m.getObjCoefficients(),
     614                         m.getRowLower(),m.getRowUpper());
     615    model.primal();
     616    int which[13] = {0,1,2,3,4,5,6,7,8,9,10,11,12};
     617    double valueIncrease[13];
     618    int sequenceIncrease[13];
     619    double valueDecrease[13];
     620    int sequenceDecrease[13];
     621    // ranging
     622    model.primalRanging(13,which,valueIncrease,sequenceIncrease,
     623                      valueDecrease,sequenceDecrease);
     624    int i;
     625    for ( i=0;i<13;i++)
     626      printf("%d increase %g %d, decrease %g %d\n",
     627             i,valueIncrease[i],sequenceIncrease[i],
     628             valueDecrease[i],sequenceDecrease[i]);
     629    assert (fabs(valueDecrease[3]-0.642857)<1.0e-4);
     630    assert (fabs(valueDecrease[8]-2.95113)<1.0e-4);
    605631  }
    606632  // test steepest edge
     
    9771003    delete [] element;
    9781004  }
    979 #if 0
     1005#if 1
    9801006  // Test barrier
    9811007  {
     
    12021228    solution.setLogLevel(63);
    12031229    solution.quadraticPrimal(1);
     1230    int numberRows = solution.numberRows();
     1231
     1232    double * rowPrimal = solution.primalRowSolution();
     1233    double * rowDual = solution.dualRowSolution();
     1234    double * rowLower = solution.rowLower();
     1235    double * rowUpper = solution.rowUpper();
     1236   
     1237    int iRow;
     1238    printf("Rows\n");
     1239    for (iRow=0;iRow<numberRows;iRow++) {
     1240      printf("%d primal %g dual %g low %g up %g\n",
     1241             iRow,rowPrimal[iRow],rowDual[iRow],
     1242             rowLower[iRow],rowUpper[iRow]);
     1243    }
     1244    double * columnPrimal = solution.primalColumnSolution();
     1245    double * columnDual = solution.dualColumnSolution();
     1246    double * columnLower = solution.columnLower();
     1247    double * columnUpper = solution.columnUpper();
    12041248    objValue = solution.getObjValue();
     1249    int iColumn;
     1250    printf("Columns\n");
     1251    for (iColumn=0;iColumn<numberColumns;iColumn++) {
     1252      printf("%d primal %g dual %g low %g up %g\n",
     1253             iColumn,columnPrimal[iColumn],columnDual[iColumn],
     1254             columnLower[iColumn],columnUpper[iColumn]);
     1255    }
    12051256    //assert(eq(objValue,3.2368421));
    12061257    //exit(77);
  • trunk/include/ClpCholeskyBase.hpp

    r328 r389  
    7979public:
    8080  virtual ClpCholeskyBase * clone() const = 0;
    81 protected:
    8281 
    8382  /// Returns type
    8483  inline int type() const
    8584  { return type_;};
     85protected:
    8686  /// Sets type
    8787  void setType(int type) {type_=type;};
     
    9393      The data members are protected to allow access for derived classes. */
    9494   //@{
    95    /// type (may be useful)
     95   /// type (may be useful) if > 20 do KKT
    9696   int type_;
    9797  /// pivotTolerance.
  • trunk/include/ClpInterior.hpp

    r388 r389  
    221221  /// Return number fixed to see if worth presolving
    222222  int numberFixed() const;
    223   /// fix variables interior says should be
    224   void fixFixed();
     223  /** fix variables interior says should be.  If reallyFix false then just
     224      set values to exact bounds */
     225  void fixFixed(bool reallyFix=true);
     226  /// Primal erturbation vector
     227  inline double * primalR() const
     228  { return primalR_;};
     229  /// Dual erturbation vector
     230  inline double * dualR() const
     231  { return dualR_;};
    225232  //@}
    226233
     
    254261  /// Checks solution
    255262  void checkSolution();
     263  /** Modifies djs to allow for quadratic.
     264      returns quadratic offset */
     265  double quadraticDjs(double * djRegion, const double * solution,
     266                      double scaleFactor);
    256267
    257268  /// To say a variable is fixed
     
    470481  double * deltaSU_;
    471482  double * deltaSL_;
     483  /// Primal regularization array
     484  double * primalR_;
     485  /// Dual regularization array
     486  double * dualR_;
    472487  /// rhs B
    473488  double * rhsB_;
  • trunk/include/ClpMatrixBase.hpp

    r372 r389  
    8787  /** Realy really scales column copy
    8888      Only called if scales already exist.
    89       Up to user ro delete */
     89      Up to user to delete */
    9090  inline virtual ClpMatrixBase * scaledColumnCopy(ClpModel * model) const
    9191  { return this->clone();};
  • trunk/include/ClpModel.hpp

    r372 r389  
    303303   inline void setRowScale(double * scale) { rowScale_ = scale;};
    304304   inline void setColumnScale(double * scale) { columnScale_ = scale;};
     305  /// Scaling of objective
     306  inline double objectiveScale() const
     307          { return objectiveScale_;} ;
     308  inline void setObjectiveScale(double value)
     309          { objectiveScale_ = value;} ;
     310  /// Scaling of rhs and bounds
     311  inline double rhsScale() const
     312          { return rhsScale_;} ;
     313  inline void setRhsScale(double value)
     314          { rhsScale_ = value;} ;
    305315   /// Sets or unsets scaling, 0 -off, 1 equilibrium, 2 geometric, 3, auto, 4 dynamic(later)
    306316   void scaling(int mode=1);
     
    560570  /// Small element value
    561571  double smallElement_;
     572  /// Scaling of objective
     573  double objectiveScale_;
     574  /// Scaling of rhs and bounds
     575  double rhsScale_;
    562576  /// Number of rows
    563577  int numberRows_;
  • trunk/include/ClpPredictorCorrector.hpp

    r373 r389  
    6161                   const double * saveRegion1, const double * saveRegion2,
    6262                   bool gentleRefine);
    63   //method: sees if looks plausible change in complementarity
    64   bool checkGoodMove(const bool doCorrector,double & bestNextGap);
     63  /// sees if looks plausible change in complementarity
     64  bool checkGoodMove(const bool doCorrector,double & bestNextGap,
     65                     bool allowIncreasingGap);
    6566  ///:  checks for one step size
    66   bool checkGoodMove2(const double move,double & bestNextGap);
     67  bool checkGoodMove2(const double move,double & bestNextGap,
     68                      bool allowIncreasingGap);
    6769  /// updateSolution.  Updates solution at end of iteration
    6870  //returns number fixed
     
    7072  ///  Save info on products of affine deltaT*deltaW and deltaS*deltaZ
    7173  double affineProduct();
     74  ///See exactly what would happen given current deltas
     75  void debugMove(int phase,double primalStep, double dualStep);
    7276  //@}
    7377
  • trunk/include/ClpQuadraticObjective.hpp

    r225 r389  
    4444                        int numberExtendedColumns_=-1);
    4545 
    46   /// Copy constructor
    47   ClpQuadraticObjective(const ClpQuadraticObjective &);
     46  /** Copy constructor .
     47      If type is -1 then make sure half symmetric,
     48      if +1 then make sure full
     49  */
     50  ClpQuadraticObjective(const ClpQuadraticObjective & rhs,int type=0);
    4851  /** Subset constructor.  Duplicates are allowed
    4952      and order is as given.
  • trunk/include/ClpSimplex.hpp

    r383 r389  
    195195                  double * costIncrease, int * sequenceIncrease,
    196196                  double * costDecrease, int * sequenceDecrease);
     197  /** Primal ranging.
     198      This computes increase/decrease in value for each given variable and corresponding
     199      sequence numbers which would change basis.  Sequence numbers are 0..numberColumns
     200      and numberColumns.. for artificials/slacks.
     201      For basic variables the sequence number will be that of the basic variables.
     202
     203      Up to user to providen correct length arrays.
     204
     205      Returns non-zero if infeasible unbounded etc
     206  */
     207  int primalRanging(int numberCheck,const int * which,
     208                  double * valueIncrease, int * sequenceIncrease,
     209                  double * valueDecrease, int * sequenceDecrease);
    197210  /// Passes in factorization
    198211  void setFactorization( ClpFactorization & factorization);
     
    525538  inline int * pivotVariable() const
    526539          { return pivotVariable_;};
    527   /// Scaling of objective
    528   inline double objectiveScale() const
    529           { return objectiveScale_;} ;
    530   inline void setObjectiveScale(double value)
    531           { objectiveScale_ = value;} ;
    532   /// Scaling of rhs and bounds
    533   inline double rhsScale() const
    534           { return rhsScale_;} ;
    535   inline void setRhsScale(double value)
    536           { rhsScale_ = value;} ;
    537540  /// If automatic scaling on
    538541  inline bool automaticScaling() const
     
    829832  /// Large bound value (for complementarity etc)
    830833  double largeValue_;
    831   /// Scaling of objective
    832   double objectiveScale_;
    833   /// Scaling of rhs and bounds
    834   double rhsScale_;
    835834  /// Largest error on Ax-b
    836835  double largestPrimalError_;
     
    10531052  /// Returns previous objective (if -1) - current if (0)
    10541053  double lastObjective(int back=1) const;
     1054  /// Set real primal infeasibility and move back
     1055  void setInfeasibility(double value);
     1056  /// Returns real primal infeasibility (if -1) - current if (0)
     1057  double lastInfeasibility(int back=1) const;
    10551058  /// Modify objective e.g. if dual infeasible in dual
    10561059  void modifyObjective(double value);
     
    10801083  /// Sum of infeasibilities for algorithm
    10811084  double infeasibility_[CLP_PROGRESS];
     1085  /// Sum of real primal infeasibilities for primal
     1086  double realInfeasibility_[CLP_PROGRESS];
    10821087#define CLP_CYCLE 12
    10831088  /// For cycle checking
  • trunk/include/ClpSimplexOther.hpp

    r343 r389  
    3939                  double * costIncrease, int * sequenceIncrease,
    4040                  double * costDecrease, int * sequenceDecrease);
     41  /** Primal ranging.
     42      This computes increase/decrease in value for each given variable and corresponding
     43      sequence numbers which would change basis.  Sequence numbers are 0..numberColumns
     44      and numberColumns.. for artificials/slacks.
     45      For basic variables the sequence number will be that of the basic variables.
     46
     47      Up to user to providen correct length arrays.
     48
     49      When here - guaranteed optimal
     50  */
     51  void primalRanging(int numberCheck,const int * which,
     52                  double * valueIncrease, int * sequenceIncrease,
     53                  double * valueDecrease, int * sequenceDecrease);
    4154  /**
    4255      Row array has row part of pivot row
     
    4457      This is used in dual ranging
    4558  */
    46   void checkRatios(CoinIndexedVector * rowArray,
     59  void checkDualRatios(CoinIndexedVector * rowArray,
    4760                   CoinIndexedVector * columnArray,
    4861                   double & costIncrease, int & sequenceIncrease,
    4962                   double & costDecrease, int & sequencedecrease);
     63  /**
     64      Row array has pivot column
     65      This is used in primal ranging
     66  */
     67  void checkPrimalRatios(CoinIndexedVector * rowArray,
     68                         int direction);
    5069  //@}
    5170};
  • trunk/include/ClpSolve.hpp

    r263 r389  
    7676      2 - interrupt handling - 0 yes, 1 no (for threadsafe)
    7777      3 - whether to make +- 1matrix - 0 yes, 1 no
     78      4 - for barrier
     79                   0 - dense cholesky
     80                   1 - Wssmp allowing some long columns
     81                   2 - Wssmp not allowing long columns
     82                   3 - Wssmp using KKT
     83                   4 - bit set to do scaling
     84                   8 - set to be aggressive with gamma/delta?
    7885  */
    7986  void setSpecialOption(int which,int value,int extraInfo=-1);
Note: See TracChangeset for help on using the changeset viewer.