Changeset 369


Ignore:
Timestamp:
May 18, 2004 10:06:56 AM (16 years ago)
Author:
forrest
Message:

improving interior point code

Location:
trunk
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • trunk/ClpCholeskyWssmp.cpp

    r328 r369  
    537537  }
    538538  bool cleanCholesky;
    539   if (model_->numberIterations()<200)
     539  if (model_->numberIterations()<2000)
    540540    cleanCholesky=true;
    541541  else
  • trunk/ClpInterior.cpp

    r335 r369  
    7777  errorRegion_(NULL),
    7878  rhsFixRegion_(NULL),
    79   updateRegion_(NULL),
    8079  upperSlack_(NULL),
    8180  lowerSlack_(NULL),
    8281  diagonal_(NULL),
    83   weights_(NULL),
    8482  solution_(NULL),
     83  workArray_(NULL),
     84  deltaX_(NULL),
     85  deltaY_(NULL),
    8586  deltaZ_(NULL),
    86   deltaW_(NULL),
    87   deltaS_(NULL),
    8887  deltaT_(NULL),
     88  deltaSU_(NULL),
     89  deltaSL_(NULL),
     90  rhsB_(NULL),
     91  rhsU_(NULL),
     92  rhsL_(NULL),
     93  rhsZ_(NULL),
     94  rhsT_(NULL),
     95  rhsC_(NULL),
    8996  zVec_(NULL),
    90   wVec_(NULL),
     97  tVec_(NULL),
    9198  cholesky_(NULL),
    9299  numberComplementarityPairs_(0),
     100  numberComplementarityItems_(0),
    93101  maximumBarrierIterations_(200),
    94102  gonePrimalFeasible_(false),
     
    154162    errorRegion_(NULL),
    155163    rhsFixRegion_(NULL),
    156     updateRegion_(NULL),
    157164    upperSlack_(NULL),
    158165    lowerSlack_(NULL),
    159166    diagonal_(NULL),
    160     weights_(NULL),
    161167    solution_(NULL),
     168    workArray_(NULL),
     169    deltaX_(NULL),
     170    deltaY_(NULL),
    162171    deltaZ_(NULL),
    163     deltaW_(NULL),
    164     deltaS_(NULL),
    165172    deltaT_(NULL),
     173    deltaSU_(NULL),
     174    deltaSL_(NULL),
     175    rhsB_(NULL),
     176    rhsU_(NULL),
     177    rhsL_(NULL),
     178    rhsZ_(NULL),
     179    rhsT_(NULL),
     180    rhsC_(NULL),
    166181    zVec_(NULL),
    167     wVec_(NULL),
     182    tVec_(NULL),
    168183    cholesky_(NULL),
    169184    numberComplementarityPairs_(0),
     185    numberComplementarityItems_(0),
    170186    maximumBarrierIterations_(200),
    171187    gonePrimalFeasible_(false),
     
    219235  errorRegion_(NULL),
    220236  rhsFixRegion_(NULL),
    221   updateRegion_(NULL),
    222237  upperSlack_(NULL),
    223238  lowerSlack_(NULL),
    224239  diagonal_(NULL),
    225   weights_(NULL),
    226240  solution_(NULL),
     241  workArray_(NULL),
     242  deltaX_(NULL),
     243  deltaY_(NULL),
    227244  deltaZ_(NULL),
    228   deltaW_(NULL),
    229   deltaS_(NULL),
    230245  deltaT_(NULL),
     246  deltaSU_(NULL),
     247  deltaSL_(NULL),
     248  rhsB_(NULL),
     249  rhsU_(NULL),
     250  rhsL_(NULL),
     251  rhsZ_(NULL),
     252  rhsT_(NULL),
     253  rhsC_(NULL),
    231254  zVec_(NULL),
    232   wVec_(NULL),
     255  tVec_(NULL),
    233256  cholesky_(NULL)
    234257{
     
    286309  errorRegion_(NULL),
    287310  rhsFixRegion_(NULL),
    288   updateRegion_(NULL),
    289311  upperSlack_(NULL),
    290312  lowerSlack_(NULL),
    291313  diagonal_(NULL),
    292   weights_(NULL),
    293314  solution_(NULL),
     315  workArray_(NULL),
     316  deltaX_(NULL),
     317  deltaY_(NULL),
    294318  deltaZ_(NULL),
    295   deltaW_(NULL),
    296   deltaS_(NULL),
    297319  deltaT_(NULL),
     320  deltaSU_(NULL),
     321  deltaSL_(NULL),
     322  rhsB_(NULL),
     323  rhsU_(NULL),
     324  rhsL_(NULL),
     325  rhsZ_(NULL),
     326  rhsT_(NULL),
     327  rhsC_(NULL),
    298328  zVec_(NULL),
    299   wVec_(NULL),
     329  tVec_(NULL),
    300330  cholesky_(NULL),
    301331  numberComplementarityPairs_(0),
     332  numberComplementarityItems_(0),
    302333  maximumBarrierIterations_(200),
    303334  gonePrimalFeasible_(false),
     
    373404  errorRegion_ = ClpCopyOfArray(rhs.errorRegion_,numberRows_);
    374405  rhsFixRegion_ = ClpCopyOfArray(rhs.rhsFixRegion_,numberRows_);
    375   updateRegion_ = ClpCopyOfArray(rhs.updateRegion_,numberRows_);
     406  deltaY_ = ClpCopyOfArray(rhs.deltaY_,numberRows_);
    376407  upperSlack_ = ClpCopyOfArray(rhs.upperSlack_,numberRows_+numberColumns_);
    377408  lowerSlack_ = ClpCopyOfArray(rhs.lowerSlack_,numberRows_+numberColumns_);
    378409  diagonal_ = ClpCopyOfArray(rhs.diagonal_,numberRows_+numberColumns_);
    379   weights_ = ClpCopyOfArray(rhs.weights_,numberRows_+numberColumns_);
     410  deltaX_ = ClpCopyOfArray(rhs.deltaX_,numberRows_+numberColumns_);
     411  deltaZ_ = ClpCopyOfArray(rhs.deltaZ_,numberRows_+numberColumns_);
     412  deltaT_ = ClpCopyOfArray(rhs.deltaT_,numberRows_+numberColumns_);
     413  deltaSU_ = ClpCopyOfArray(rhs.deltaSU_,numberRows_+numberColumns_);
     414  deltaSL_ = ClpCopyOfArray(rhs.deltaSL_,numberRows_+numberColumns_);
     415  rhsB_ = ClpCopyOfArray(rhs.rhsB_,numberRows_);
     416  rhsU_ = ClpCopyOfArray(rhs.rhsU_,numberRows_+numberColumns_);
     417  rhsL_ = ClpCopyOfArray(rhs.rhsL_,numberRows_+numberColumns_);
     418  rhsZ_ = ClpCopyOfArray(rhs.rhsZ_,numberRows_+numberColumns_);
     419  rhsT_ = ClpCopyOfArray(rhs.rhsT_,numberRows_+numberColumns_);
     420  rhsC_ = ClpCopyOfArray(rhs.rhsC_,numberRows_+numberColumns_);
    380421  solution_ = ClpCopyOfArray(rhs.solution_,numberRows_+numberColumns_);
    381   deltaZ_ = ClpCopyOfArray(rhs.deltaZ_,numberRows_+numberColumns_);
    382   deltaW_ = ClpCopyOfArray(rhs.deltaW_,numberRows_+numberColumns_);
    383   deltaS_ = ClpCopyOfArray(rhs.deltaS_,numberRows_+numberColumns_);
    384   deltaT_ = ClpCopyOfArray(rhs.deltaT_,numberRows_+numberColumns_);
     422  workArray_ = ClpCopyOfArray(rhs.workArray_,numberRows_+numberColumns_);
    385423  zVec_ = ClpCopyOfArray(rhs.zVec_,numberRows_+numberColumns_);
    386   wVec_ = ClpCopyOfArray(rhs.wVec_,numberRows_+numberColumns_);
     424  tVec_ = ClpCopyOfArray(rhs.tVec_,numberRows_+numberColumns_);
    387425  cholesky_ = rhs.cholesky_->clone();
    388426  numberComplementarityPairs_ = rhs.numberComplementarityPairs_;
     427  numberComplementarityItems_ = rhs.numberComplementarityItems_;
    389428  maximumBarrierIterations_ = rhs.maximumBarrierIterations_;
    390429  gonePrimalFeasible_ = rhs.gonePrimalFeasible_;
     
    424463  delete [] rhsFixRegion_;
    425464  rhsFixRegion_ = NULL;
    426   delete [] updateRegion_;
    427   updateRegion_ = NULL;
     465  delete [] deltaY_;
     466  deltaY_ = NULL;
    428467  delete [] upperSlack_;
    429468  upperSlack_ = NULL;
     
    432471  delete [] diagonal_;
    433472  diagonal_ = NULL;
    434   delete [] weights_;
    435   weights_ = NULL;
     473  delete [] deltaX_;
     474  deltaX_ = NULL;
     475  delete [] deltaZ_;
     476  deltaZ_ = NULL;
     477  delete [] deltaT_;
     478  deltaT_ = NULL;
     479  delete [] deltaSU_;
     480  deltaSU_ = NULL;
     481  delete [] deltaSL_;
     482  deltaSL_ = NULL;
     483  delete [] rhsB_;
     484  rhsB_ = NULL;
     485  delete [] rhsU_;
     486  rhsU_ = NULL;
     487  delete [] rhsL_;
     488  rhsL_ = NULL;
     489  delete [] rhsZ_;
     490  rhsZ_ = NULL;
     491  delete [] rhsT_;
     492  rhsT_ = NULL;
     493  delete [] rhsC_;
     494  rhsC_ = NULL;
    436495  delete [] solution_;
    437496  solution_ = NULL;
    438   delete [] deltaZ_;
    439   deltaZ_ = NULL;
    440   delete [] deltaW_;
    441   deltaW_ = NULL;
    442   delete [] deltaS_;
    443   deltaS_ = NULL;
    444   delete [] deltaT_;
    445   deltaT_ = NULL;
     497  delete [] workArray_;
     498  workArray_ = NULL;
    446499  delete [] zVec_;
    447500  zVec_ = NULL;
    448   delete [] wVec_;
    449   wVec_ = NULL;
     501  delete [] tVec_;
     502  tVec_ = NULL;
    450503  delete cholesky_;
    451504}
     
    510563  assert (!rhsFixRegion_);
    511564  rhsFixRegion_ = new double [numberRows_];
    512   assert (!updateRegion_);
    513   updateRegion_ = new double [numberRows_];
     565  assert (!deltaY_);
     566  deltaY_ = new double [numberRows_];
     567  CoinZeroN(deltaY_,numberRows_);
    514568  assert (!upperSlack_);
    515569  upperSlack_ = new double [nTotal];
     
    518572  assert (!diagonal_);
    519573  diagonal_ = new double [nTotal];
    520   assert (!weights_);
    521   weights_ = new double [nTotal];
     574  assert (!deltaX_);
     575  deltaX_ = new double [nTotal];
     576  CoinZeroN(deltaX_,nTotal);
    522577  assert (!deltaZ_);
    523578  deltaZ_ = new double [nTotal];
    524579  CoinZeroN(deltaZ_,nTotal);
    525   assert (!deltaW_);
    526   deltaW_ = new double [nTotal];
    527   CoinZeroN(deltaW_,nTotal);
    528   assert (!deltaS_);
    529   deltaS_ = new double [nTotal];
    530580  assert (!deltaT_);
    531581  deltaT_ = new double [nTotal];
     582  CoinZeroN(deltaT_,nTotal);
     583  assert (!deltaSU_);
     584  deltaSU_ = new double [nTotal];
     585  CoinZeroN(deltaSU_,nTotal);
     586  assert (!deltaSL_);
     587  deltaSL_ = new double [nTotal];
     588  CoinZeroN(deltaSL_,nTotal);
     589  assert (!rhsB_);
     590  rhsB_ = new double [numberRows_];
     591  CoinZeroN(rhsB_,numberRows_);
     592  assert (!rhsU_);
     593  rhsU_ = new double [nTotal];
     594  CoinZeroN(rhsU_,nTotal);
     595  assert (!rhsL_);
     596  rhsL_ = new double [nTotal];
     597  CoinZeroN(rhsL_,nTotal);
     598  assert (!rhsZ_);
     599  rhsZ_ = new double [nTotal];
     600  CoinZeroN(rhsZ_,nTotal);
     601  assert (!rhsT_);
     602  rhsT_ = new double [nTotal];
     603  CoinZeroN(rhsT_,nTotal);
     604  assert (!rhsC_);
     605  rhsC_ = new double [nTotal];
     606  CoinZeroN(rhsC_,nTotal);
     607  assert (!workArray_);
     608  workArray_ = new double [nTotal];
     609  CoinZeroN(workArray_,nTotal);
    532610  assert (!zVec_);
    533611  zVec_ = new double [nTotal];
    534612  CoinZeroN(zVec_,nTotal);
    535   assert (!wVec_);
    536   wVec_ = new double [nTotal];
    537   CoinZeroN(wVec_,nTotal);
     613  assert (!tVec_);
     614  tVec_ = new double [nTotal];
     615  CoinZeroN(tVec_,nTotal);
    538616  assert (!dj_);
    539617  dj_ = new double [nTotal];
     
    566644  delete [] rhsFixRegion_;
    567645  rhsFixRegion_ = NULL;
    568   delete [] updateRegion_;
    569   updateRegion_ = NULL;
     646  delete [] deltaY_;
     647  deltaY_ = NULL;
    570648  delete [] upperSlack_;
    571649  upperSlack_ = NULL;
     
    574652  delete [] diagonal_;
    575653  diagonal_ = NULL;
    576   delete [] weights_;
    577   weights_ = NULL;
    578   delete [] deltaZ_;
    579   deltaZ_ = NULL;
    580   delete [] deltaW_;
    581   deltaW_ = NULL;
    582   delete [] deltaS_;
    583   deltaS_ = NULL;
    584   delete [] deltaT_;
    585   deltaT_ = NULL;
     654  delete [] deltaX_;
     655  deltaX_ = NULL;
     656  delete [] workArray_;
     657  workArray_ = NULL;
    586658  delete [] zVec_;
    587659  zVec_ = NULL;
    588   delete [] wVec_;
    589   wVec_ = NULL;
     660  delete [] tVec_;
     661  tVec_ = NULL;
    590662  delete [] dj_;
    591663  dj_ = NULL;
  • trunk/ClpPredictorCorrector.cpp

    r365 r369  
    99
    1010#include "CoinPragma.hpp"
    11 static int xxxxxx=-1;
    1211#include <math.h>
    1312
     
    6261  //set iterations
    6362  numberIterations_=-1;
     63  int numberTotal = numberRows_+numberColumns_;
    6464  //initialize solution here
    6565  if(createSolution()<0) {
     
    9898  double * savePi=NULL;
    9999  double * savePrimal=NULL;
    100   int numberTotal = numberRows_+numberColumns_;
    101100  // Extra regions for centering
    102   double * saveWeights = new double[numberTotal];
     101  double * saveX = new double[numberTotal];
     102  double * saveY = new double[numberRows_];
    103103  double * saveZ = new double[numberTotal];
    104   double * saveW = new double[numberTotal];
     104  double * saveT = new double[numberTotal];
    105105  while (problemStatus_<0) {
    106     if (numberIterations_==xxxxxx) {
    107       FILE *fp = fopen("yy.yy","wb");
    108       if (fp) {
    109         int nrow=numberRows_;
    110         int ncol=numberColumns_;
    111         // and flip
    112         int i;
    113         for (i=0;i<nrow;i++) {
    114           solution_[ncol+i]= -solution_[ncol+i];
    115           double value;
    116           value=lowerSlack_[ncol+i];
    117           lowerSlack_[ncol+i]=upperSlack_[ncol+i];
    118           upperSlack_[ncol+i]=value;
    119           value=zVec_[ncol+i];
    120           zVec_[ncol+i]=wVec_[ncol+i];
    121           wVec_[ncol+i]=value;
    122         }
    123         fwrite(dual_,sizeof(double),nrow,fp);
    124         fwrite(errorRegion_,sizeof(double),nrow,fp);
    125         fwrite(rhsFixRegion_,sizeof(double),nrow,fp);
    126         fwrite(solution_+ncol,sizeof(double),nrow,fp);
    127         fwrite(diagonal_+ncol,sizeof(double),nrow,fp);
    128         fwrite(wVec_+ncol,sizeof(double),nrow,fp);
    129         fwrite(zVec_+ncol,sizeof(double),nrow,fp);
    130         fwrite(upperSlack_+ncol,sizeof(double),nrow,fp);
    131         fwrite(lowerSlack_+ncol,sizeof(double),nrow,fp);
    132         fwrite(solution_,sizeof(double),ncol,fp);
    133         fwrite(diagonal_,sizeof(double),ncol,fp);
    134         fwrite(wVec_,sizeof(double),ncol,fp);
    135         fwrite(zVec_,sizeof(double),ncol,fp);
    136         fwrite(upperSlack_,sizeof(double),ncol,fp);
    137         fwrite(lowerSlack_,sizeof(double),ncol,fp);
    138         fclose(fp);
    139         actualDualStep_=0.0;
    140         actualPrimalStep_=0.0;
    141       } else {
    142         printf("no file\n");
    143       }
    144     }
    145106    //#define FULL_DEBUG
    146107#ifdef FULL_DEBUG
     
    153114      printf(" col  dsol  ddiag  dwvec  dzvec dbdslu dbdsll\n");
    154115      for (i=0;i<numberColumns_+numberRows_;i++) {
    155         printf(" %d %g %g %g %g %g %g\n",i,solution_[i],diagonal_[i],wVec_[i],
     116        printf(" %d %g %g %g %g %g %g\n",i,solution_[i],diagonal_[i],tVec_[i],
    156117               zVec_[i],upperSlack_[i],lowerSlack_[i]);
    157118      }
    158119    }
    159120#endif
    160     complementarityGap_=complementarityGap(numberComplementarityPairs_,0);
     121    complementarityGap_=complementarityGap(numberComplementarityPairs_,
     122                                           numberComplementarityItems_,0);
    161123      handler_->message(CLP_BARRIER_ITERATION,messages_)
    162124    <<numberIterations_
     
    311273      }
    312274    }
    313     bool goodMove=false;
    314     double bestNextGap=COIN_DBL_MAX;
    315     double nextGap=bestNextGap;
     275    double nextGap=COIN_DBL_MAX;
    316276    int nextNumber=0;
     277    int nextNumberItems=0;
    317278    worstDirectionAccuracy_=0.0;
    318     while (!goodMove) {
    319       goodMove=true;
    320       int newDropped=0;
    321       //Predictor step
    322       //prepare for cholesky.  Set up scaled diagonal in weights
    323       //  ** for efficiency may be better if scale factor known before
    324       double norm2=0.0;
    325       double maximumValue;
    326       getNorms(diagonal_,numberTotal,maximumValue,norm2);
    327       diagonalNorm_ = sqrt(norm2/numberComplementarityPairs_);
    328       diagonalScaleFactor_=1.0;
    329       double maximumAllowable=eScale;
    330       //scale so largest is less than allowable ? could do better
    331       double factor=0.5;
    332       while (maximumValue>maximumAllowable) {
    333         diagonalScaleFactor_*=factor;
    334         maximumValue*=factor;
    335       } /* endwhile */
    336       if (diagonalScaleFactor_!=1.0) {
    337         handler_->message(CLP_BARRIER_SCALING,messages_)
    338           <<"diagonal"<<diagonalScaleFactor_
    339           <<CoinMessageEol;
    340         diagonalNorm_*=diagonalScaleFactor_;
    341       }
    342       multiplyAdd(NULL,numberTotal,0.0,diagonal_,
    343                   diagonalScaleFactor_);
    344       int * rowsDroppedThisTime = new int [numberRows_];
    345       newDropped=cholesky_->factorize(diagonal_,rowsDroppedThisTime);
    346       if (newDropped>0) {
    347         //int newDropped2=cholesky_->factorize(diagonal_,rowsDroppedThisTime);
    348         //assert(!newDropped2);
    349         if (newDropped<0) {
    350           //replace dropped
    351           newDropped=-newDropped;
    352           //off 1 to allow for reset all
    353           newDropped--;
    354           //set all bits false
    355           cholesky_->resetRowsDropped();
    356         }
    357       } else if (newDropped==-1) {
    358         printf("Out of memory\n");
    359         problemStatus_=4;
    360         return -1;
    361       }
    362       delete [] rowsDroppedThisTime;
    363       if (cholesky_->status()) {
    364         std::cout << "bad cholesky?" <<std::endl;
    365         abort();
    366       }
    367       int phase=0; // predictor, corrector , primal dual
    368       double directionAccuracy=0.0;
    369       bool doCorrector=true;
    370       //goodMove=false;
    371       if (goodMove) {
    372         //set up for affine direction
    373         setupForSolve(phase);
    374         directionAccuracy=findDirectionVector(phase);
    375         if (directionAccuracy>worstDirectionAccuracy_) {
    376           worstDirectionAccuracy_=directionAccuracy;
    377         }
    378         findStepLength(phase,NULL);
    379         nextGap=complementarityGap(nextNumber,1);
    380        
    381         bestNextGap=max(nextGap,0.8*complementarityGap_);
    382         if (complementarityGap_>1.0e-5*numberComplementarityPairs_) {
    383           //std::cout <<"predicted duality gap "<<nextGap<<std::endl;
    384           double part1=nextGap/numberComplementarityPairs_;
    385           double part2=nextGap/complementarityGap_;
    386           mu_=part1*part2*part2;
     279    int newDropped=0;
     280    //Predictor step
     281    //prepare for cholesky.  Set up scaled diagonal in deltaX
     282    //  ** for efficiency may be better if scale factor known before
     283    double norm2=0.0;
     284    double maximumValue;
     285    getNorms(diagonal_,numberTotal,maximumValue,norm2);
     286    diagonalNorm_ = sqrt(norm2/numberComplementarityPairs_);
     287    diagonalScaleFactor_=1.0;
     288    double maximumAllowable=eScale;
     289    //scale so largest is less than allowable ? could do better
     290    double factor=0.5;
     291    while (maximumValue>maximumAllowable) {
     292      diagonalScaleFactor_*=factor;
     293      maximumValue*=factor;
     294    } /* endwhile */
     295    if (diagonalScaleFactor_!=1.0) {
     296      handler_->message(CLP_BARRIER_SCALING,messages_)
     297        <<"diagonal"<<diagonalScaleFactor_
     298        <<CoinMessageEol;
     299      diagonalNorm_*=diagonalScaleFactor_;
     300    }
     301    multiplyAdd(NULL,numberTotal,0.0,diagonal_,
     302                diagonalScaleFactor_);
     303    int * rowsDroppedThisTime = new int [numberRows_];
     304    newDropped=cholesky_->factorize(diagonal_,rowsDroppedThisTime);
     305    if (newDropped>0) {
     306      //int newDropped2=cholesky_->factorize(diagonal_,rowsDroppedThisTime);
     307      //assert(!newDropped2);
     308      if (newDropped<0) {
     309        //replace dropped
     310        newDropped=-newDropped;
     311        //off 1 to allow for reset all
     312        newDropped--;
     313        //set all bits false
     314        cholesky_->resetRowsDropped();
     315      }
     316    } else if (newDropped==-1) {
     317      printf("Out of memory\n");
     318      problemStatus_=4;
     319      return -1;
     320    }
     321    delete [] rowsDroppedThisTime;
     322    if (cholesky_->status()) {
     323      std::cout << "bad cholesky?" <<std::endl;
     324      abort();
     325    }
     326    int phase=0; // predictor, corrector , primal dual
     327    double directionAccuracy=0.0;
     328    bool doCorrector=true;
     329    bool goodMove=true;
     330    //set up for affine direction
     331    setupForSolve(phase);
     332    directionAccuracy=findDirectionVector(phase);
     333    if (directionAccuracy>worstDirectionAccuracy_) {
     334      worstDirectionAccuracy_=directionAccuracy;
     335    }
     336    findStepLength(phase);
     337    nextGap=complementarityGap(nextNumber,nextNumberItems,1);
     338    double affineGap=nextGap;
     339    int bestPhase=0;
     340    double bestNextGap=nextGap;
     341    // ?
     342    bestNextGap=max(nextGap,0.8*complementarityGap_);
     343    if (complementarityGap_>1.0e-4*numberComplementarityPairs_) {
     344      //std::cout <<"predicted duality gap "<<nextGap<<std::endl;
     345      double part1=nextGap/numberComplementarityPairs_;
     346      double part2=nextGap/complementarityGap_;
     347      mu_=part1*part2*part2;
    387348#if 0
    388           double papermu =complementarityGap_/numberComplementarityPairs_;
    389           double affmu = nextGap/nextNumber;
    390           double sigma = pow(affmu/papermu,3);
    391           printf("mu %g, papermu %g, affmu %g, sigma %g sigmamu %g\n",
    392                 mu_,papermu,affmu,sigma,sigma*papermu);
     349      double papermu =complementarityGap_/numberComplementarityPairs_;
     350      double affmu = nextGap/nextNumber;
     351      double sigma = pow(affmu/papermu,3);
     352      printf("mu %g, papermu %g, affmu %g, sigma %g sigmamu %g\n",
     353            mu_,papermu,affmu,sigma,sigma*papermu);
    393354#endif 
    394           //printf("paper mu %g\n",(nextGap*nextGap*nextGap)/(complementarityGap_*complementarityGap_*
    395           //                                        (double) numberComplementarityPairs_));
    396         } else {
    397           double phi;
    398           if (numberComplementarityPairs_<=500) {
    399             phi=pow((double) numberComplementarityPairs_,2.0);
    400           } else {
    401             phi=pow((double) numberComplementarityPairs_,1.5);
    402             if (phi<500.0*500.0) {
    403               phi=500.0*500.0;
    404             }
    405           }
    406           mu_=complementarityGap_/phi;
    407           //? could use gapO?
    408           //if small then should be stopping
    409           //if (mu_<1.0e-4/(numberComplementarityPairs_*qqqq)) {
    410           //mu_=1.0e-4/(numberComplementarityPairs_*qqqq);
    411           //? better to skip corrector?
    412           //}
    413         }
    414         //save information
    415         double product=affineProduct();
    416         //can we do corrector step?
    417         double xx= complementarityGap_*(beta2-tau) +product;
    418         //std::cout<<"gap part "<<
    419         //complementarityGap_*(beta2-tau)<<" after adding product = "<<xx<<std::endl;
    420         //xx=-1.0;
    421         if (xx>0.0) {
    422           double saveMu = mu_;
    423           double mu2=numberComplementarityPairs_;
    424           mu2=xx/mu2;
    425           if (mu2>mu_) {
    426             //std::cout<<" could increase to "<<mu2<<std::endl;
    427             //was mu2=mu2*0.25;
    428             mu2=mu2*0.99;
    429             if (mu2<mu_) {
    430               mu_=mu2;
    431               //std::cout<<" changing mu to "<<mu_<<std::endl;
    432             } else {
    433               //std::cout<<std::endl;
    434             }
    435           } else {
    436             //std::cout<<" should decrease to "<<mu2<<std::endl;
    437             mu_=0.5*mu2;
    438             //std::cout<<" changing mu to "<<mu_<<std::endl;
    439           }
    440           handler_->message(CLP_BARRIER_MU,messages_)
    441             <<saveMu<<mu_
    442             <<CoinMessageEol;
    443         } else {
    444           //std::cout<<" bad by any standards"<<std::endl;
    445         }
    446         //printf("product %g mu %g\n",product,mu_);
    447         if (complementarityGap_*(beta2-tau)+product-mu_*numberComplementarityPairs_<0.0) {
    448           doCorrector=false;
    449           goodMove=false;
    450           bestNextGap=COIN_DBL_MAX;
    451           handler_->message(CLP_BARRIER_INFO,messages_)
    452             <<"no corrector step"
    453             <<CoinMessageEol;
    454         } else {
    455           phase=1;
     355      //printf("paper mu %g\n",(nextGap*nextGap*nextGap)/(complementarityGap_*complementarityGap_*
     356      //                                            (double) numberComplementarityPairs_));
     357    } else {
     358      double phi;
     359      if (numberComplementarityPairs_<=500) {
     360        phi=pow((double) numberComplementarityPairs_,2.0);
     361      } else {
     362        phi=pow((double) numberComplementarityPairs_,1.5);
     363        if (phi<500.0*500.0) {
     364          phi=500.0*500.0;
    456365        }
    457366      }
    458       if (goodMove&&doCorrector) {
    459         //set up for next step
    460         setupForSolve(phase);
    461         double directionAccuracy2=findDirectionVector(phase);
    462         if (directionAccuracy2>worstDirectionAccuracy_) {
    463           worstDirectionAccuracy_=directionAccuracy2;
     367      mu_=complementarityGap_/phi;
     368    }
     369    //save information
     370    double product=affineProduct();
     371    //#define ALWAYS
     372#ifndef ALWAYS
     373    //can we do corrector step?
     374    double xx= complementarityGap_*(beta2-tau) +product;
     375    if (xx>0.0) {
     376      double saveMu = mu_;
     377      double mu2=numberComplementarityPairs_;
     378      mu2=xx/mu2;
     379      if (mu2>mu_) {
     380        //std::cout<<" could increase to "<<mu2<<std::endl;
     381        //was mu2=mu2*0.25;
     382        mu2=mu2*0.99;
     383        if (mu2<mu_) {
     384          mu_=mu2;
     385          //std::cout<<" changing mu to "<<mu_<<std::endl;
     386        } else {
     387          //std::cout<<std::endl;
    464388        }
    465         double testValue=1.0e2*directionAccuracy;
    466         if (1.0e2*projectionTolerance_>testValue) {
    467           testValue=1.0e2*projectionTolerance_;
     389      } else {
     390        //std::cout<<" should decrease to "<<mu2<<std::endl;
     391        mu_=0.5*mu2;
     392        //std::cout<<" changing mu to "<<mu_<<std::endl;
     393      }
     394      handler_->message(CLP_BARRIER_MU,messages_)
     395        <<saveMu<<mu_
     396        <<CoinMessageEol;
     397    } else {
     398      //std::cout<<" bad by any standards"<<std::endl;
     399    }
     400    if (complementarityGap_*(beta2-tau)+product-mu_*numberComplementarityPairs_<0.0) {
     401#ifdef SOME_DEBUG
     402      printf("failed 1 product %g mu %g - %g < 0.0, nextGap %g\n",product,mu_,
     403             complementarityGap_*(beta2-tau)+product-mu_*numberComplementarityPairs_,
     404             nextGap);
     405#endif
     406      doCorrector=false;
     407      if (nextGap>0.9*complementarityGap_||1) {
     408        goodMove=false;
     409        bestNextGap=COIN_DBL_MAX;
     410      }
     411      //double floatNumber = 2.0*numberComplementarityPairs_;
     412      //floatNumber = 1.0*numberComplementarityItems_;
     413      //mu_=nextGap/floatNumber;
     414      handler_->message(CLP_BARRIER_INFO,messages_)
     415        <<"no corrector step"
     416        <<CoinMessageEol;
     417    } else {
     418      phase=1;
     419    }
     420#else
     421    phase=1;
     422#endif
     423    if (goodMove&&doCorrector) {
     424      //set up for next step
     425      setupForSolve(phase);
     426      double directionAccuracy2=findDirectionVector(phase);
     427      if (directionAccuracy2>worstDirectionAccuracy_) {
     428        worstDirectionAccuracy_=directionAccuracy2;
     429      }
     430      double testValue=1.0e2*directionAccuracy;
     431      if (1.0e2*projectionTolerance_>testValue) {
     432        testValue=1.0e2*projectionTolerance_;
     433      }
     434      if (primalTolerance()>testValue) {
     435        testValue=primalTolerance();
     436      }
     437      if (maximumRHSError_>testValue) {
     438        testValue=maximumRHSError_;
     439      }
     440      if (directionAccuracy2>testValue&&numberIterations_>=-77) {
     441        goodMove=false;
     442#ifdef SOME_DEBUG
     443        printf("accuracy %g phase 1 failed, test value %g\n",
     444               directionAccuracy2,testValue);
     445#endif
     446        bestNextGap=COIN_DBL_MAX;
     447      }
     448      if (goodMove) {
     449        findStepLength(phase);
     450        // just for debug
     451        nextGap = complementarityGap(nextNumber,nextNumberItems,1);
     452#ifndef ALWAYS
     453        if (numberIterations_>=-77) {
     454          goodMove=checkGoodMove(true,bestNextGap);
     455          if (!goodMove) {
     456#ifdef SOME_DEBUG
     457            printf("checkGoodMove failed\n");
     458#endif
     459            if (affineGap<0.5*complementarityGap_) {
     460              //||(affineGap<0.95*complementarityGap_&&(numberIterations_%3)==0)) {
     461              // Back to affine
     462              setupForSolve(0);
     463              directionAccuracy=findDirectionVector(0);
     464              findStepLength(0);
     465              nextGap=complementarityGap(nextNumber,nextNumberItems,1);
     466              goodMove=true;
     467            }
     468          }
     469        } else {
     470          goodMove=true;
    468471        }
    469         if (primalTolerance()>testValue) {
    470           testValue=primalTolerance();
    471         }
    472         if (maximumRHSError_>testValue) {
    473           testValue=maximumRHSError_;
    474         }
    475         if (directionAccuracy2>testValue&&numberIterations_>=-77) {
    476           goodMove=false;
    477           bestNextGap=COIN_DBL_MAX;
    478         }
    479         if (goodMove) {
    480           findStepLength(phase,NULL);
    481           // just for debug
    482           nextGap = complementarityGap(nextNumber,1);
    483           if (numberIterations_>=-77) {
    484             goodMove=checkGoodMove(true,bestNextGap);
    485           } else {
    486             goodMove=true;
    487           }
    488         }
    489       }
    490       if (!goodMove) {
    491         // Just primal dual step
    492         goodMove=true;
    493         double floatNumber;
    494         floatNumber = 2.0*numberComplementarityPairs_;
    495         //floatNumber = 1.1*numberComplementarityPairs_;
    496         mu_=complementarityGap_/floatNumber;
    497         //set up for next step
    498         setupForSolve(2);
    499         double directionAccuracy=findDirectionVector(2);
    500         if (directionAccuracy>worstDirectionAccuracy_) {
    501           worstDirectionAccuracy_=directionAccuracy;
    502         }
    503         double testValue=1.0e2*directionAccuracy;
    504         if (1.0e2*projectionTolerance_>testValue) {
    505           testValue=1.0e2*projectionTolerance_;
    506         }
    507         if (primalTolerance()>testValue) {
    508           testValue=primalTolerance();
    509         }
    510         if (maximumRHSError_>testValue) {
    511           testValue=maximumRHSError_;
    512         }
    513         findStepLength(2,NULL);
    514         // just for debug
    515         nextGap=complementarityGap(nextNumber,2);
     472#endif
    516473      }
    517     } /* endwhile */
    518     //numberFixed=updateSolution();
    519     //numberFixedTotal+=numberFixed;
    520     if (numberIterations_==1) {
    521       FILE *fp = fopen("xx.xx","rb");
    522       if (fp) {
    523         int nrow=numberRows_;
    524         int ncol=numberColumns_;
    525         fread(dual_,sizeof(double),nrow,fp);
    526         fread(errorRegion_,sizeof(double),nrow,fp);
    527         fread(rhsFixRegion_,sizeof(double),nrow,fp);
    528         fread(solution_+ncol,sizeof(double),nrow,fp);
    529         fread(diagonal_+ncol,sizeof(double),nrow,fp);
    530         fread(wVec_+ncol,sizeof(double),nrow,fp);
    531         fread(zVec_+ncol,sizeof(double),nrow,fp);
    532         fread(upperSlack_+ncol,sizeof(double),nrow,fp);
    533         fread(lowerSlack_+ncol,sizeof(double),nrow,fp);
    534         fread(solution_,sizeof(double),ncol,fp);
    535         fread(diagonal_,sizeof(double),ncol,fp);
    536         fread(wVec_,sizeof(double),ncol,fp);
    537         fread(zVec_,sizeof(double),ncol,fp);
    538         fread(upperSlack_,sizeof(double),ncol,fp);
    539         fread(lowerSlack_,sizeof(double),ncol,fp);
    540         fclose(fp);
    541         // and flip
    542         int i;
    543         for (i=0;i<nrow;i++) {
    544           solution_[ncol+i]= -solution_[ncol+i];
    545           double value;
    546           value=lowerSlack_[ncol+i];
    547           lowerSlack_[ncol+i]=upperSlack_[ncol+i];
    548           upperSlack_[ncol+i]=value;
    549           value=zVec_[ncol+i];
    550           zVec_[ncol+i]=wVec_[ncol+i];
    551           wVec_[ncol+i]=value;
    552         }
    553         actualDualStep_=0.0;
    554         actualPrimalStep_=0.0;
    555       } else {
    556         printf("no file\n");
     474    }
     475    if (!goodMove) {
     476      // Just primal dual step
     477      double floatNumber;
     478      floatNumber = 2.0*numberComplementarityPairs_;
     479      //floatNumber = 1.1*numberComplementarityPairs_;
     480      mu_=complementarityGap_/floatNumber;
     481      //set up for next step
     482      setupForSolve(2);
     483      double directionAccuracy=findDirectionVector(2);
     484      if (directionAccuracy>worstDirectionAccuracy_) {
     485        worstDirectionAccuracy_=directionAccuracy;
     486      }
     487      double testValue=1.0e2*directionAccuracy;
     488      if (1.0e2*projectionTolerance_>testValue) {
     489        testValue=1.0e2*projectionTolerance_;
     490      }
     491      if (primalTolerance()>testValue) {
     492        testValue=primalTolerance();
     493      }
     494      if (maximumRHSError_>testValue) {
     495        testValue=maximumRHSError_;
     496      }
     497      findStepLength(2);
     498      // just for debug
     499      nextGap=complementarityGap(nextNumber,nextNumberItems,2);
     500      if (nextGap>bestNextGap&&bestPhase==0&&0) {
     501        // Back to affine
     502        setupForSolve(0);
     503        directionAccuracy=findDirectionVector(0);
     504        findStepLength(0);
     505        nextGap=complementarityGap(nextNumber,nextNumberItems,1);
    557506      }
    558507    }
    559 #if 0
    560     // See if we can do centering steps
    561     memcpy(saveWeights,weights_,numberTotal*sizeof(double));
    562     memcpy(saveZ,deltaZ_,numberTotal*sizeof(double));
    563     memcpy(saveW,deltaW_,numberTotal*sizeof(double));
    564     double savePrimalStep = actualPrimalStep_;
    565     double saveDualStep = actualDualStep_;
    566     double saveMu = mu_;
    567     mu_=nextGap / (double) nextNumber;
    568     setupForSolve(3);
    569     findDirectionVector(3);
    570     memcpy(deltaZ_,saveZ,numberTotal*sizeof(double));
    571     memcpy(deltaW_,saveW,numberTotal*sizeof(double));
    572     findStepLength(3,saveWeights);
    573     double xGap = complementarityGap(nextNumber,3);
    574 #if 0
    575     goodMove=checkGoodMove(true,bestNextGap);
    576 #endif
    577     if (xGap>nextGap) {
    578       mu_=saveMu;
    579       actualPrimalStep_ = savePrimalStep;
    580       actualDualStep_ = saveDualStep;
    581       memcpy(weights_,saveWeights,numberTotal*sizeof(double));
    582       memcpy(deltaZ_,saveZ,numberTotal*sizeof(double));
    583       memcpy(deltaW_,saveW,numberTotal*sizeof(double));
     508    if (!goodMove)
     509      mu_=nextGap / ((double) 1.1*nextNumber);
     510    goodMove=true;
     511    // Do centering steps
     512    int numberTries=0;
     513    double nextCenterGap=0.0;
     514    int numberGoodTries=0;
     515    double originalDualStep=actualDualStep_;
     516    double originalPrimalStep=actualPrimalStep_;
     517    while (goodMove&&numberTries<4) {
     518      goodMove=false;
     519      numberTries++;
     520      memcpy(saveX,deltaX_,numberTotal*sizeof(double));
     521      memcpy(saveY,deltaY_,numberRows_*sizeof(double));
     522      memcpy(saveZ,deltaZ_,numberTotal*sizeof(double));
     523      memcpy(saveT,deltaT_,numberTotal*sizeof(double));
     524      double savePrimalStep = actualPrimalStep_;
     525      double saveDualStep = actualDualStep_;
     526      double saveMu = mu_;
     527      setupForSolve(3);
     528      findDirectionVector(3);
     529      findStepLength(3);
     530      double xGap = complementarityGap(nextNumber,nextNumberItems,3);
     531      // If one small then that's the one that counts
     532      double checkDual=saveDualStep;
     533      double checkPrimal=savePrimalStep;
     534      if (checkDual>5.0*checkPrimal) {
     535        checkDual=2.0*checkPrimal;
     536      } else if (checkPrimal>5.0*checkDual) {
     537        checkPrimal=2.0*checkDual;
     538      }
     539      //if (actualPrimalStep_<checkPrimal||
     540      //  actualDualStep_<checkDual||
     541      //  (xGap>nextGap&&xGap>0.9*complementarityGap_)) {
     542      if (actualPrimalStep_<checkPrimal||
     543          actualDualStep_<checkDual) {
     544#ifdef SOME_DEBUG
     545        printf("PP rejected gap %g, steps %g %g, 2 gap %g, steps %g %g\n",xGap,
     546               actualPrimalStep_,actualDualStep_,nextGap,savePrimalStep,saveDualStep);
     547#endif
     548        mu_=saveMu;
     549        actualPrimalStep_ = savePrimalStep;
     550        actualDualStep_ = saveDualStep;
     551        memcpy(deltaX_,saveX,numberTotal*sizeof(double));
     552        memcpy(deltaY_,saveY,numberRows_*sizeof(double));
     553        memcpy(deltaZ_,saveZ,numberTotal*sizeof(double));
     554        memcpy(deltaT_,saveT,numberTotal*sizeof(double));
     555      } else {
     556#ifdef SOME_DEBUG
     557        printf("PPphase 3 gap %g, steps %g %g, 2 gap %g, steps %g %g\n",xGap,
     558               actualPrimalStep_,actualDualStep_,nextGap,savePrimalStep,saveDualStep);
     559#endif
     560        numberGoodTries++;
     561        nextCenterGap=xGap;
     562        // See if big enough change
     563        if (actualPrimalStep_<1.01*checkPrimal||
     564            actualDualStep_<1.01*checkDual) {
     565          // stop now
     566        } else {
     567          // carry on
     568          goodMove=true;
     569        }
     570      }
    584571    }
    585 #endif
     572    if (numberGoodTries&&handler_->logLevel()>1) {
     573      printf("%d centering steps moved from (gap %g, dual %g, primal %g) to (gap %g, dual %g, primal %g)\n",
     574             numberGoodTries,nextGap,originalDualStep,originalPrimalStep,
     575             nextCenterGap, actualDualStep_,actualPrimalStep_);
     576    }
    586577    numberFixed=updateSolution();
    587578    numberFixedTotal+=numberFixed;
    588579  } /* endwhile */
    589   delete [] saveWeights;
     580  delete [] saveX;
     581  delete [] saveY;
    590582  delete [] saveZ;
    591   delete [] saveW;
     583  delete [] saveT;
    592584  if (savePi) {
    593585    //std::cout<<"Restoring from iteration "<<saveIteration<<std::endl;
     
    615607    <<objectiveValue()
    616608    <<CoinMessageEol;
     609  //#ifdef SOME_DEBUG
     610  if (handler_->logLevel()>1)
     611    printf("ENDRUN status %d after %d iterations\n",problemStatus_,numberIterations_);
     612  //#endif
    617613  //std::cout<<"Absolute primal infeasibility at end "<<sumPrimalInfeasibilities_<<std::endl;
    618614  //std::cout<<"Absolute dual infeasibility at end "<<sumDualInfeasibilities_<<std::endl;
     
    628624//         1 corrector
    629625//         2 primal dual
    630 double ClpPredictorCorrector::findStepLength(const int phase,const double * oldWeight)
     626double ClpPredictorCorrector::findStepLength( int phase)
    631627{
    632628  double directionNorm=0.0;
     
    637633  int chosenPrimalSequence=-1;
    638634  int chosenDualSequence=-1;
    639   double extra=eExtra;
    640635  double * zVec = zVec_;
    641   double * wVec = wVec_;
     636  double * tVec = tVec_;
    642637  double * primal = solution_;
    643   double * lower = lower_;
    644   double * upper = upper_;
    645638  double * lowerSlack = lowerSlack_;
    646639  double * upperSlack = upperSlack_;
    647   double * work1 = deltaZ_;
    648   double * work2 = deltaW_;
    649   double * work3 = deltaS_;
    650   double * work4 = deltaT_;
    651   //direction vector in weights
    652   double * weights = weights_;
     640  //direction vector in deltaX
     641  double * deltaX = deltaX_;
    653642  int iColumn;
    654   switch (phase) {
    655   case 0:
    656     //Now get affine deltas for Z(duals on LBds) and W (duals on UBds)
    657     for (iColumn=0;iColumn<numberTotal;iColumn++) {
    658       if (!flagged(iColumn)) {
    659         double z1=-zVec[iColumn];
    660         double w1=-wVec[iColumn];
    661         double work3Value=0.0;
    662         double work4Value=0.0;
    663         double directionElement=weights[iColumn];
    664         double value=primal[iColumn];
    665         if (directionNorm<fabs(directionElement)) {
    666           directionNorm=fabs(directionElement);
    667         }
    668         //below does not feel right - can it be simplified because
    669         //of zero values for zVec and wVec
    670         if (lowerBound(iColumn)||
    671                         upperBound(iColumn)) {
    672           if (lowerBound(iColumn)) {
    673             double gap=lowerSlack[iColumn]+extra;
    674             double delta;
    675             if (!fakeLower(iColumn)) {
    676               delta = lower[iColumn]+lowerSlack[iColumn]
    677              -  value - directionElement;
    678             } else {
    679               delta = - directionElement;
    680             }
    681             z1+=(zVec[iColumn]*delta)/gap;
    682             work3Value=-delta;
    683             if (lowerSlack[iColumn]<maximumPrimalStep*delta) {
    684               maximumPrimalStep=lowerSlack[iColumn]/delta;
    685               chosenPrimalSequence=iColumn;
    686             }
    687             if (zVec[iColumn]>tolerance) {
    688               if (zVec[iColumn]<-z1*maximumDualStep) {
    689                 maximumDualStep=-zVec[iColumn]/z1;
    690                 chosenDualSequence=iColumn;
    691               }
    692             }
    693           }
    694           if (upperBound(iColumn)) {
    695             double gap=upperSlack[iColumn]+extra;
    696             double delta;
    697             if (!fakeUpper(iColumn)) {
    698               delta = -upper[iColumn]+upperSlack[iColumn]
    699                         + value + directionElement;
    700             } else {
    701               delta = directionElement;
    702             }
    703             w1+=(wVec[iColumn]*delta)/gap;
    704             work4Value=-delta;
    705             if (upperSlack[iColumn]<maximumPrimalStep*delta) {
    706               maximumPrimalStep=upperSlack[iColumn]/delta;
    707               chosenPrimalSequence=iColumn;
    708             }
    709             if (wVec[iColumn]>tolerance) {
    710               if (wVec[iColumn]<-w1*maximumDualStep) {
    711                 maximumDualStep=-wVec[iColumn]/w1;
    712                 chosenDualSequence=iColumn;
    713               }
    714             }
    715           }
    716         } else {
    717           //free
    718           double gap=fabs(value);
    719           double multiplier=1.0/gap;
    720           if (gap<1.0) {
    721             multiplier=1,0;
    722           }
    723           z1-=multiplier*directionElement*zVec[iColumn];
    724           w1+=multiplier*directionElement*wVec[iColumn];
    725         }
    726         work1[iColumn]=z1;
    727         work2[iColumn]=w1;
    728         work3[iColumn]=work3Value;
    729         work4[iColumn]=work4Value;
     643  for (iColumn=0;iColumn<numberTotal;iColumn++) {
     644    if (!flagged(iColumn)) {
     645      double z1=0.0;
     646      double t1=0.0;
     647      double directionElement=deltaX[iColumn];
     648      double value=primal[iColumn];
     649      if (directionNorm<fabs(directionElement)) {
     650        directionNorm=fabs(directionElement);
     651      }
     652      if (lowerBound(iColumn)||
     653          upperBound(iColumn)) {
     654        if (lowerBound(iColumn)) {
     655          double delta = - deltaSL_[iColumn];
     656          z1 = deltaZ_[iColumn];
     657          if (lowerSlack[iColumn]<maximumPrimalStep*delta) {
     658            maximumPrimalStep=lowerSlack[iColumn]/delta;
     659            chosenPrimalSequence=iColumn;
     660          }
     661          if (zVec[iColumn]>tolerance) {
     662            if (zVec[iColumn]<-z1*maximumDualStep) {
     663              maximumDualStep=-zVec[iColumn]/z1;
     664              chosenDualSequence=iColumn;
     665            }
     666          }
     667        } else {
     668          assert (!zVec[iColumn]);
     669        }
     670        if (upperBound(iColumn)) {
     671          double delta = - deltaSU_[iColumn];;
     672          t1 = deltaT_[iColumn];
     673          if (upperSlack[iColumn]<maximumPrimalStep*delta) {
     674            maximumPrimalStep=upperSlack[iColumn]/delta;
     675            chosenPrimalSequence=iColumn;
     676          }
     677          if (tVec[iColumn]>tolerance) {
     678            if (tVec[iColumn]<-t1*maximumDualStep) {
     679              maximumDualStep=-tVec[iColumn]/t1;
     680              chosenDualSequence=iColumn;
     681            }
     682          }
     683        } else {
     684          assert (!tVec[iColumn]);
     685        }
    730686      } else {
    731         work1[iColumn]=0.0;
    732         work2[iColumn]=0.0;
    733         work3[iColumn]=0.0;
    734         work4[iColumn]=0.0;
    735       }
    736     }
    737     break;
    738   case 1:
    739     //corrector step
    740     for ( iColumn=0;iColumn<numberTotal;iColumn++) {
    741       if (!flagged(iColumn)) {
    742         double z1=-zVec[iColumn];
    743         double w1=-wVec[iColumn];
    744         double work3Value=0.0;
    745         double work4Value=0.0;
    746         double directionElement=weights[iColumn];
    747         double value=primal[iColumn];
    748         if (directionNorm<fabs(directionElement)) {
    749           directionNorm=fabs(directionElement);
    750         }
    751         //below does not feel right - can it be simplified because
    752         //of zero values for zVec and wVec
    753         if (lowerBound(iColumn)||
    754                         upperBound(iColumn)) {
    755           if (lowerBound(iColumn)) {
    756             double gap=lowerSlack[iColumn]+extra;
    757             double delta;
    758             if (!fakeLower(iColumn)) {
    759               delta = lower[iColumn]+lowerSlack[iColumn]
    760              -  value - directionElement;
    761             } else {
    762               delta = - directionElement;
    763             }
    764             z1+=(mu_-work3[iColumn]+zVec[iColumn]*delta)/gap;
    765             work3Value=-delta;
    766             if (lowerSlack[iColumn]<maximumPrimalStep*delta) {
    767               maximumPrimalStep=lowerSlack[iColumn]/delta;
    768               chosenPrimalSequence=iColumn;
    769             }
    770             if (zVec[iColumn]>tolerance) {
    771               if (zVec[iColumn]<-z1*maximumDualStep) {
    772                 maximumDualStep=-zVec[iColumn]/z1;
    773                 chosenDualSequence=iColumn;
    774               }
    775             }
    776           }
    777           if (upperBound(iColumn)) {
    778             double gap=upperSlack[iColumn]+extra;
    779             double delta;
    780             if (!fakeUpper(iColumn)) {
    781               delta = -upper[iColumn]+upperSlack[iColumn]
    782                         + value + directionElement;
    783             } else {
    784               delta = directionElement;
    785             }
    786             //double delta = -upper[iColumn]+upperSlack[iColumn]-
    787               //+ value + directionElement;
    788             w1+=(mu_-work4[iColumn]+wVec[iColumn]*delta)/gap;
    789             work4Value=-delta;
    790             if (upperSlack[iColumn]<maximumPrimalStep*delta) {
    791               maximumPrimalStep=upperSlack[iColumn]/delta;
    792               chosenPrimalSequence=iColumn;
    793             }
    794             if (wVec[iColumn]>tolerance) {
    795               if (wVec[iColumn]<-w1*maximumDualStep) {
    796                 maximumDualStep=-wVec[iColumn]/w1;
    797                 chosenDualSequence=iColumn;
    798               }
    799             }
    800           }
    801         } else {
    802           //free
    803           double gap=fabs(value);
    804           double multiplier=1.0/gap;
    805           if (gap<1.0) {
    806             multiplier=1,0;
    807           }
    808           z1+=multiplier*(mu_-work3[iColumn]-directionElement*zVec[iColumn]);
    809           w1+=multiplier*(mu_-work4[iColumn]+directionElement*wVec[iColumn]);
    810         }
    811         work1[iColumn]=z1;
    812         work2[iColumn]=w1;
    813         work3[iColumn]=work3Value;
    814         work4[iColumn]=work4Value;
    815       } else {
    816         work1[iColumn]=0.0;
    817         work2[iColumn]=0.0;
    818         work3[iColumn]=0.0;
    819         work4[iColumn]=0.0;
    820       }
    821     }
    822     break;
    823   case 2:
    824     //just primal dual
    825     for ( iColumn=0;iColumn<numberTotal;iColumn++) {
    826       if (!flagged(iColumn)) {
    827         double z1=-zVec[iColumn];
    828         double w1=-wVec[iColumn];
    829         double work3Value=0.0;
    830         double work4Value=0.0;
    831         double directionElement=weights[iColumn];
    832         double value=primal[iColumn];
    833         if (directionNorm<fabs(directionElement)) {
    834           directionNorm=fabs(directionElement);
    835         }
    836         //below does not feel right - can it be simplified because
    837         //of zero values for zVec and wVec
    838         if (lowerBound(iColumn)||
    839                         upperBound(iColumn)) {
    840           if (lowerBound(iColumn)) {
    841             double gap=lowerSlack[iColumn]+extra;
    842             double delta;
    843             if (!fakeLower(iColumn)) {
    844               delta = lower[iColumn]+lowerSlack[iColumn]
    845              -  value - directionElement;
    846             } else {
    847               delta = - directionElement;
    848             }
    849             z1+=(mu_+zVec[iColumn]*delta)/gap;
    850             work3Value=-delta;
    851             if (lowerSlack[iColumn]<maximumPrimalStep*delta) {
    852               maximumPrimalStep=lowerSlack[iColumn]/delta;
    853               chosenPrimalSequence=iColumn;
    854             }
    855             if (zVec[iColumn]>tolerance) {
    856               if (zVec[iColumn]<-z1*maximumDualStep) {
    857                 maximumDualStep=-zVec[iColumn]/z1;
    858                 chosenDualSequence=iColumn;
    859               }
    860             }
    861           }
    862           if (upperBound(iColumn)) {
    863             double gap=upperSlack[iColumn]+extra;
    864             double delta;
    865             if (!fakeUpper(iColumn)) {
    866               delta = -upper[iColumn]+upperSlack[iColumn]
    867                         + value + directionElement;
    868             } else {
    869               delta = directionElement;
    870             }
    871             w1+=(mu_+wVec[iColumn]*delta)/gap;
    872             work4Value=-delta;
    873             if (upperSlack[iColumn]<maximumPrimalStep*delta) {
    874               maximumPrimalStep=upperSlack[iColumn]/delta;
    875               chosenPrimalSequence=iColumn;
    876             }
    877             if (wVec[iColumn]>tolerance) {
    878               if (wVec[iColumn]<-w1*maximumDualStep) {
    879                 maximumDualStep=-wVec[iColumn]/w1;
    880                 chosenDualSequence=iColumn;
    881               }
    882             }
    883           }
    884         } else {
    885           //free
    886           double gap=fabs(value);
    887           double multiplier=1.0/gap;
    888           if (gap<1.0) {
    889             multiplier=1,0;
    890           }
    891           z1+=multiplier*(mu_-directionElement*zVec[iColumn]);
    892           w1+=multiplier*(mu_+directionElement*wVec[iColumn]);
    893         }
    894         work1[iColumn]=z1;
    895         work2[iColumn]=w1;
    896         work3[iColumn]=work3Value;
    897         work4[iColumn]=work4Value;
    898       } else {
    899         work1[iColumn]=0.0;
    900         work2[iColumn]=0.0;
    901         work3[iColumn]=0.0;
    902         work4[iColumn]=0.0;
    903       }
    904     }
    905     break;
    906   case 3:
    907     //centering step
    908     for ( iColumn=0;iColumn<numberTotal;iColumn++) {
    909       if (!flagged(iColumn)) {
    910         // old deltas from previous step
    911         double z1=work1[iColumn];
    912         double w1=work2[iColumn];
    913         double work3Value=0.0;
    914         double work4Value=0.0;
    915         double directionElement=weights[iColumn];
    916         double oldElement = oldWeight[iColumn];
    917         double newElement = oldElement+directionElement;
    918         weights[iColumn]=newElement;
    919         double value=primal[iColumn];
    920         if (directionNorm<fabs(directionElement)) {
    921           directionNorm=fabs(directionElement);
    922         }
    923         if (lowerBound(iColumn)||
    924                         upperBound(iColumn)) {
    925           if (lowerBound(iColumn)) {
    926             double gap=lowerSlack[iColumn]+extra;
    927             double delta;
    928             double zDelta;
    929             if (!fakeLower(iColumn)) {
    930               zDelta = lower[iColumn]+lowerSlack[iColumn]
    931              -  value - directionElement;
    932             } else {
    933               zDelta = - directionElement;
    934             }
    935             delta = zDelta - oldElement;
    936             z1+=(work3[iColumn]-zVec[iColumn]*zDelta)/gap;
    937             work3Value=-delta;
    938             if (lowerSlack[iColumn]<maximumPrimalStep*delta) {
    939               maximumPrimalStep=lowerSlack[iColumn]/delta;
    940               chosenPrimalSequence=iColumn;
    941             }
    942             if (zVec[iColumn]>tolerance) {
    943               if (zVec[iColumn]<-z1*maximumDualStep) {
    944                 maximumDualStep=-zVec[iColumn]/z1;
    945                 chosenDualSequence=iColumn;
    946               }
    947             }
    948           }
    949           if (upperBound(iColumn)) {
    950             double gap=upperSlack[iColumn]+extra;
    951             double delta;
    952             double zDelta;
    953             if (!fakeUpper(iColumn)) {
    954               zDelta = -upper[iColumn]+upperSlack[iColumn]
    955                         + value + directionElement;
    956             } else {
    957               zDelta = directionElement;
    958             }
    959             delta = zDelta + oldElement;
    960             //double delta = -upper[iColumn]+upperSlack[iColumn]-
    961               //+ value + directionElement;
    962             w1+=(work4[iColumn]+wVec[iColumn]*zDelta)/gap;
    963             work4Value=-delta;
    964             if (upperSlack[iColumn]<maximumPrimalStep*delta) {
    965               maximumPrimalStep=upperSlack[iColumn]/delta;
    966               chosenPrimalSequence=iColumn;
    967             }
    968             if (wVec[iColumn]>tolerance) {
    969               if (wVec[iColumn]<-w1*maximumDualStep) {
    970                 maximumDualStep=-wVec[iColumn]/w1;
    971                 chosenDualSequence=iColumn;
    972               }
    973             }
    974           }
    975         } else {
    976           //free
    977           double gap=fabs(value);
    978           double multiplier=1.0/gap;
    979           if (gap<1.0) {
    980             multiplier=1,0;
    981           }
    982           z1+=multiplier*(work3[iColumn]-directionElement*zVec[iColumn]);
    983           w1+=multiplier*(work4[iColumn]+directionElement*wVec[iColumn]);
    984         }
    985         work1[iColumn]=z1;
    986         work2[iColumn]=w1;
    987         work3[iColumn]=work3Value;
    988         work4[iColumn]=work4Value;
    989       } else {
    990         work1[iColumn]=0.0;
    991         work2[iColumn]=0.0;
    992         work3[iColumn]=0.0;
    993         work4[iColumn]=0.0;
    994       }
    995     }
    996     break;
    997   }
    998   //printf("step - phase %d, norm %g, dual step %g, primal step %g\n",
    999   // phase,directionNorm,maximumPrimalStep,maximumDualStep);
     687        //free
     688        double gap=fabs(value);
     689        double multiplier=1.0/gap;
     690        if (gap<1.0) {
     691          multiplier=1,0;
     692        }
     693        z1 = -zVec[iColumn]-multiplier*directionElement*zVec[iColumn];
     694        t1 = -tVec[iColumn] + multiplier*directionElement*tVec[iColumn];
     695        deltaZ_[iColumn]=z1;
     696        deltaT_[iColumn]=t1;
     697      }
     698    }
     699  }
     700#ifdef SOME_DEBUG
     701  printf("new step - phase %d, norm %g, dual step %g, primal step %g\n",
     702         phase,directionNorm,maximumPrimalStep,maximumDualStep);
     703#endif
    1000704  actualPrimalStep_=stepLength_*maximumPrimalStep;
    1001705  if (phase>=0&&actualPrimalStep_>1.0) {
     
    1006710    actualDualStep_=1.0;
    1007711  }
     712#ifdef FULL_DEBUG
     713  if (phase==3){
     714    double minBeta = 0.1*mu_;
     715    double maxBeta = 10.0*mu_;
     716    for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
     717      if (!flagged(iColumn)) {
     718        if (lowerBound(iColumn)) {
     719          double change = rhsL_[iColumn] + deltaX_[iColumn];
     720          double dualValue=zVec[iColumn]+actualDualStep_*deltaZ_[iColumn];
     721          double primalValue=lowerSlack[iColumn]+actualPrimalStep_*change;
     722          double gapProduct=dualValue*primalValue;
     723          if (delta2Z_[iColumn]<minBeta||delta2Z_[iColumn]>maxBeta)
     724            printf("3lower %d primal %g, dual %g, gap %g, old gap %g\n",
     725                   iColumn,primalValue,dualValue,gapProduct,delta2Z_[iColumn]);
     726        } 
     727        if (upperBound(iColumn)) {
     728          double change = rhsU_[iColumn]-deltaX_[iColumn];
     729          double dualValue=tVec[iColumn]+actualDualStep_*deltaT_[iColumn];
     730          double primalValue=upperSlack[iColumn]+actualPrimalStep_*change;
     731          double gapProduct=dualValue*primalValue;
     732          if (delta2T_[iColumn]<minBeta||delta2T_[iColumn]>maxBeta)
     733            printf("3upper %d primal %g, dual %g, gap %g, old gap %g\n",
     734                 iColumn,primalValue,dualValue,gapProduct,delta2T_[iColumn]);
     735        }
     736      }
     737    }
     738  }
     739#endif
     740#ifdef SOME_DEBUG
     741  {
     742    double largestL=0.0;
     743    double smallestL=COIN_DBL_MAX;
     744    double largestU=0.0;
     745    double smallestU=COIN_DBL_MAX;
     746    double sumL=0.0;
     747    double sumU=0.0;
     748    int nL=0;
     749    int nU=0;
     750    for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
     751      if (!flagged(iColumn)) {
     752        if (lowerBound(iColumn)) {
     753          double change = rhsL_[iColumn] + deltaX_[iColumn];
     754          double dualValue=zVec[iColumn]+actualDualStep_*deltaZ_[iColumn];
     755          double primalValue=lowerSlack[iColumn]+actualPrimalStep_*change;
     756          double gapProduct=dualValue*primalValue;
     757          largestL = max(largestL,gapProduct);
     758          smallestL = min(smallestL,gapProduct);
     759          nL++;
     760          sumL += gapProduct;
     761        } 
     762        if (upperBound(iColumn)) {
     763          double change = rhsU_[iColumn]-deltaX_[iColumn];
     764          double dualValue=tVec[iColumn]+actualDualStep_*deltaT_[iColumn];
     765          double primalValue=upperSlack[iColumn]+actualPrimalStep_*change;
     766          double gapProduct=dualValue*primalValue;
     767          largestU = max(largestU,gapProduct);
     768          smallestU = min(smallestU,gapProduct);
     769          nU++;
     770          sumU += gapProduct;
     771        }
     772      }
     773    }
     774    double mu = (sumL+sumU)/((double) (nL+nU));
     775
     776    double minBeta = 0.1*mu;
     777    double maxBeta = 10.0*mu;
     778    int nBL=0;
     779    int nAL=0;
     780    int nBU=0;
     781    int nAU=0;
     782    for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
     783      if (!flagged(iColumn)) {
     784        if (lowerBound(iColumn)) {
     785          double change = rhsL_[iColumn] + deltaX_[iColumn];
     786          double dualValue=zVec[iColumn]+actualDualStep_*deltaZ_[iColumn];
     787          double primalValue=lowerSlack[iColumn]+actualPrimalStep_*change;
     788          double gapProduct=dualValue*primalValue;
     789          if (gapProduct<minBeta)
     790            nBL++;
     791          else if (gapProduct>maxBeta)
     792            nAL++;
     793        } 
     794        if (upperBound(iColumn)) {
     795          double change = rhsU_[iColumn]-deltaX_[iColumn];
     796          double dualValue=tVec[iColumn]+actualDualStep_*deltaT_[iColumn];
     797          double primalValue=upperSlack[iColumn]+actualPrimalStep_*change;
     798          double gapProduct=dualValue*primalValue;
     799          if (gapProduct<minBeta)
     800            nBU++;
     801          else if (gapProduct>maxBeta)
     802            nAU++;
     803        }
     804      }
     805    }
     806    printf("phase %d new mu %g new gap %g\n",phase,mu,sumL+sumU);
     807    printf("          %d lower, smallest %g, %d below - largest %g, %d above\n",
     808           nL,smallestL,nBL,largestL,nAL);
     809    printf("          %d upper, smallest %g, %d below - largest %g, %d above\n",
     810           nU,smallestU,nBU,largestU,nAU);
     811  }
     812#endif
    1008813  return directionNorm;
    1009814}
     
    1081886double ClpPredictorCorrector::findDirectionVector(const int phase)
    1082887{
    1083   if (phase==3) {
    1084     // centering step
    1085     double * work2 = deltaW_;
    1086     int numberTotal = numberRows_+numberColumns_;
    1087     double * tempRegion = new double [numberRows_];
    1088     double * newError = new double [numberRows_];
    1089     // For KKT separate out
    1090 #ifndef KKT
    1091     int iColumn;
    1092     for (iColumn=0;iColumn<numberTotal;iColumn++)
    1093       weights_[iColumn] = work2[iColumn];
    1094     multiplyAdd(weights_+numberColumns_,numberRows_,-1.0,tempRegion,0.0);
    1095     matrix_->times(1.0,weights_,tempRegion);
    1096 #else
    1097     // regions in will be work2 and newError
    1098     // regions out weights_ and tempRegion
    1099     memset(newError,0,numberRows_*sizeof(double));
    1100     double * region1Save=NULL;//for refinement
    1101 #endif
    1102 #ifndef KKT
    1103     double maximumRHS = maximumAbsElement(tempRegion,numberRows_);
    1104     double scale=1.0;
    1105     double unscale=1.0;
    1106     if (maximumRHS>1.0e-30) {
    1107       if (maximumRHS<=0.5) {
    1108         double factor=2.0;
    1109         while (maximumRHS<=0.5) {
    1110           maximumRHS*=factor;
    1111           scale*=factor;
    1112         } /* endwhile */
    1113       } else if (maximumRHS>=2.0&&maximumRHS<=COIN_DBL_MAX) {
    1114         double factor=0.5;
    1115         while (maximumRHS>=2.0) {
    1116           maximumRHS*=factor;
    1117           scale*=factor;
    1118         } /* endwhile */
    1119       }
    1120       unscale=diagonalScaleFactor_/scale;
    1121     } else {
    1122       //effectively zero
    1123       scale=0.0;
    1124       unscale=0.0;
    1125     }
    1126     //printf("--putting scales to 1.0\n");
    1127     //scale=1.0;
    1128     //unscale=1.0;
    1129     multiplyAdd(NULL,numberRows_,0.0,tempRegion,scale);
    1130     cholesky_->solve(tempRegion);
    1131     multiplyAdd(NULL,numberRows_,0.0,tempRegion,unscale);
    1132     //CoinZeroN(newError,numberRows_);
    1133     multiplyAdd(tempRegion,numberRows_,-1.0,weights_+numberColumns_,0.0);
    1134     CoinZeroN(weights_,numberColumns_);
    1135     matrix_->transposeTimes(1.0,tempRegion,weights_);
    1136     //if flagged then entries zero so can do
    1137     for (iColumn=0;iColumn<numberTotal;iColumn++)
    1138       weights_[iColumn] = weights_[iColumn]*diagonal_[iColumn]
    1139         -work2[iColumn];
    1140 #else
    1141     solveSystem(weights_, tempRegion,
    1142                 work2,newError,region1Save,regionSave,lastError>1.0e-5);
    1143 #endif
    1144     multiplyAdd(weights_+numberColumns_,numberRows_,-1.0,newError,0.0);
    1145     matrix_->times(1.0,weights_,newError);
    1146     delete [] newError;
    1147     delete [] tempRegion;
    1148     return 0.0;
    1149   }
    1150888  double projectionTolerance=projectionTolerance_;
    1151889  //temporary
     
    1162900  }
    1163901  double * newError = new double [numberRows_];
    1164   double * work2 = deltaW_;
     902  double * workArray = workArray_;
    1165903  int numberTotal = numberRows_+numberColumns_;
    1166904  //if flagged then entries zero so can do
     
    1169907  int iColumn;
    1170908  for (iColumn=0;iColumn<numberTotal;iColumn++)
    1171     weights_[iColumn] = work2[iColumn] - solution_[iColumn];
    1172   multiplyAdd(weights_+numberColumns_,numberRows_,-1.0,updateRegion_,0.0);
    1173   matrix_->times(1.0,weights_,updateRegion_);
     909    deltaX_[iColumn] = workArray[iColumn] - solution_[iColumn];
     910  multiplyAdd(deltaX_+numberColumns_,numberRows_,-1.0,deltaY_,0.0);
     911  matrix_->times(1.0,deltaX_,deltaY_);
    1174912#else
    1175   // regions in will be work2 and newError
    1176   // regions out weights_ and updateRegion_
     913  // regions in will be workArray and newError
     914  // regions out deltaX_ and deltaY_
    1177915  multiplyAdd(solution_+numberColumns_,numberRows_,1.0,newError,0.0);
    1178916  matrix_->times(-1.0,solution_,newError);
     
    1188926    goodSolve=true;
    1189927#ifndef KKT
    1190     double maximumRHS = maximumAbsElement(updateRegion_,numberRows_);
     928    double maximumRHS = maximumAbsElement(deltaY_,numberRows_);
    1191929    double saveMaximum = maximumRHS;
    1192930    double scale=1.0;
     
    1215953    //scale=1.0;
    1216954    //unscale=1.0;
    1217     multiplyAdd(NULL,numberRows_,0.0,updateRegion_,scale);
    1218     cholesky_->solve(updateRegion_);
    1219     multiplyAdd(NULL,numberRows_,0.0,updateRegion_,unscale);
     955    multiplyAdd(NULL,numberRows_,0.0,deltaY_,scale);
     956    cholesky_->solve(deltaY_);
     957    multiplyAdd(NULL,numberRows_,0.0,deltaY_,unscale);
    1220958    if (numberTries) {
    1221959      //refine?
     
    1223961      if (lastError>1.0e-5)
    1224962        scaleX=0.8;
    1225       multiplyAdd(regionSave,numberRows_,1.0,updateRegion_,scaleX);
     963      multiplyAdd(regionSave,numberRows_,1.0,deltaY_,scaleX);
    1226964    }
    1227965    //CoinZeroN(newError,numberRows_);
    1228     multiplyAdd(updateRegion_,numberRows_,-1.0,weights_+numberColumns_,0.0);
    1229     CoinZeroN(weights_,numberColumns_);
    1230     matrix_->transposeTimes(1.0,updateRegion_,weights_);
     966    multiplyAdd(deltaY_,numberRows_,-1.0,deltaX_+numberColumns_,0.0);
     967    CoinZeroN(deltaX_,numberColumns_);
     968    matrix_->transposeTimes(1.0,deltaY_,deltaX_);
    1231969    //if flagged then entries zero so can do
    1232970    for (iColumn=0;iColumn<numberTotal;iColumn++)
    1233       weights_[iColumn] = weights_[iColumn]*diagonal_[iColumn]
    1234         -work2[iColumn];
     971      deltaX_[iColumn] = deltaX_[iColumn]*diagonal_[iColumn]
     972        -workArray[iColumn];
    1235973#else
    1236     solveSystem(weights_, updateRegion_,
    1237                 work2,newError,region1Save,regionSave,lastError>1.0e-5);
    1238 #endif
    1239     multiplyAdd(weights_+numberColumns_,numberRows_,-1.0,newError,0.0);
    1240     matrix_->times(1.0,weights_,newError);
     974    solveSystem(deltaX_, deltaY_,
     975                workArray,newError,region1Save,regionSave,lastError>1.0e-5);
     976#endif
     977    multiplyAdd(deltaX_+numberColumns_,numberRows_,-1.0,newError,0.0);
     978    matrix_->times(1.0,deltaX_,newError);
    1241979    numberTries++;
    1242980   
     
    12691007        newError[iRow]=result;
    12701008        //newError[iRow]=0.0;
    1271         //assert(updateRegion_[iRow]==0.0);
    1272         updateRegion_[iRow]=0.0;
     1009        //assert(deltaY_[iRow]==0.0);
     1010        deltaY_[iRow]=0.0;
    12731011      }
    12741012    }
     
    12971035#endif
    12981036        }
    1299         CoinMemcpyN(updateRegion_,numberRows_,regionSave);
     1037        CoinMemcpyN(deltaY_,numberRows_,regionSave);
    13001038#ifndef KKT
    1301         multiplyAdd(newError,numberRows_,-1.0,updateRegion_,0.0);
     1039        multiplyAdd(newError,numberRows_,-1.0,deltaY_,0.0);
    13021040#else
    1303         CoinMemcpyN(weights_,numberTotal,region1Save);
     1041        CoinMemcpyN(deltaX_,numberTotal,region1Save);
    13041042        // and back to input region
    1305         CoinMemcpyN(updateRegion_,numberRows_,newError);
     1043        CoinMemcpyN(deltaY_,numberRows_,newError);
    13061044#endif
    13071045      }
     
    13101048      //bring back previous
    13111049      relativeError=lastError;
    1312       CoinMemcpyN(regionSave,numberRows_,updateRegion_);
     1050      CoinMemcpyN(regionSave,numberRows_,deltaY_);
    13131051#ifndef KKT
    1314       multiplyAdd(updateRegion_,numberRows_,-1.0,weights_+numberColumns_,0.0);
    1315       CoinZeroN(weights_,numberColumns_);
    1316       matrix_->transposeTimes(1.0,updateRegion_,weights_);
     1052      multiplyAdd(deltaY_,numberRows_,-1.0,deltaX_+numberColumns_,0.0);
     1053      CoinZeroN(deltaX_,numberColumns_);
     1054      matrix_->transposeTimes(1.0,deltaY_,deltaX_);
    13171055      //if flagged then entries zero so can do
    13181056      for (iColumn=0;iColumn<numberTotal;iColumn++)
    1319         weights_[iColumn] = weights_[iColumn]*diagonal_[iColumn]
    1320           -work2[iColumn];
     1057        deltaX_[iColumn] = deltaX_[iColumn]*diagonal_[iColumn]
     1058          -workArray[iColumn];
    13211059#else
    1322         CoinMemcpyN(region1Save,numberTotal,weights_);
     1060        CoinMemcpyN(region1Save,numberTotal,deltaX_);
    13231061#endif
    13241062    }
     
    13291067#endif
    13301068  delete [] newError;
     1069  // now rest
     1070  double * zVec = zVec_;
     1071  double * tVec = tVec_;
     1072  double * lowerSlack = lowerSlack_;
     1073  double * upperSlack = upperSlack_;
     1074  double extra=eExtra;
     1075
     1076  for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
     1077    deltaSU_[iColumn]=0.0;
     1078    deltaSL_[iColumn]=0.0;
     1079    deltaZ_[iColumn]=0.0;
     1080    deltaT_[iColumn]=0.0;
     1081    if (!flagged(iColumn)) {
     1082      double deltaX = deltaX_[iColumn];
     1083      if (lowerBound(iColumn)||upperBound(iColumn)) {
     1084        if (lowerBound(iColumn)) {
     1085          double deltaSL = rhsL_[iColumn]+deltaX;
     1086          double slack = lowerSlack[iColumn]+extra;
     1087          deltaSL_[iColumn]=deltaSL;
     1088          deltaZ_[iColumn]=(rhsZ_[iColumn]-zVec[iColumn]*deltaSL)/slack;
     1089        }
     1090        if (upperBound(iColumn)) {
     1091          double deltaSU = rhsU_[iColumn]-deltaX;
     1092          double slack = upperSlack[iColumn]+extra;
     1093          deltaSU_[iColumn]=deltaSU;
     1094          deltaT_[iColumn]=(rhsT_[iColumn]-tVec[iColumn]*deltaSU)/slack;
     1095        }
     1096      } else {
     1097        // free
     1098        double gap=fabs(solution_[iColumn]);
     1099        double multiplier=1.0/gap;
     1100        if (gap<1.0) {
     1101          multiplier=1,0;
     1102        }
     1103        deltaZ_[iColumn] = -zVec[iColumn] - multiplier*deltaX*zVec[iColumn];
     1104        deltaT_[iColumn] = -tVec[iColumn] + multiplier*deltaX*tVec[iColumn];
     1105      }
     1106    }
     1107  }
    13311108  return relativeError;
    13321109}
     
    13731150  double infiniteCheck=1.0e40;
    13741151  //double     fakeCheck=1.0e10;
    1375   //use weights region for work region
     1152  //use deltaX region for work region
    13761153  for (iColumn=0;iColumn<numberTotal;iColumn++) {
    13771154    double primalValue = solution_[iColumn];
     
    13851162      dj_[iColumn]=0.0;
    13861163      diagonal_[iColumn]=1.0;
    1387       weights_[iColumn]=1.0;
     1164      deltaX_[iColumn]=1.0;
    13881165      double lowerValue=lower_[iColumn];
    13891166      double upperValue=upper_[iColumn];
     
    14431220      solution_[iColumn]=lower_[iColumn];
    14441221      diagonal_[iColumn]=0.0;
    1445       weights_[iColumn]=0.0;
     1222      deltaX_[iColumn]=0.0;
    14461223    }
    14471224  }
     
    14761253  cholesky_->solve(errorRegion_);
    14771254  //create information for solution
    1478   multiplyAdd(errorRegion_,numberRows_,-1.0,weights_+numberColumns_,0.0);
    1479   CoinZeroN(weights_,numberColumns_);
    1480   matrix_->transposeTimes(1.0,errorRegion_,weights_);
     1255  multiplyAdd(errorRegion_,numberRows_,-1.0,deltaX_+numberColumns_,0.0);
     1256  CoinZeroN(deltaX_,numberColumns_);
     1257  matrix_->transposeTimes(1.0,errorRegion_,deltaX_);
    14811258#else
    14821259  // reverse sign on solution
    14831260  multiplyAdd(NULL,numberRows_+numberColumns_,0.0,solution_,-1.0);
    1484   solveSystem(weights_,errorRegion_,solution_,NULL,NULL,NULL,false);
     1261  solveSystem(deltaX_,errorRegion_,solution_,NULL,NULL,NULL,false);
    14851262#endif
    14861263  //do reduced costs
     
    14961273    dj_[iColumn]=value;
    14971274  }
    1498   initialValue=1.0e2;
     1275  initialValue=1.0e1;
    14991276  if (rhsNorm_*1.0e-2>initialValue) {
    15001277    initialValue=rhsNorm_*1.0e-2;
    15011278  }
    15021279  double smallestBoundDifference=COIN_DBL_MAX;
    1503   double * fakeSolution = weights_;
     1280  double * fakeSolution = deltaX_;
    15041281  for ( iColumn=0;iColumn<numberTotal;iColumn++) {
    15051282    if (!flagged(iColumn)) {
     
    15231300  double largeGap=1.0e15;
    15241301  double safeObjectiveValue=2.0*objectiveNorm_;
    1525   //safeObjectiveValue=1.0+objectiveNorm_;
    1526   //printf("temp safe obj value of %g\n",safeObjectiveValue);
    15271302  double safeFree=1.0e-1*initialValue;
    1528   safeFree=1.0;
     1303  //safeFree=1.0;
     1304  //printf("normal safe dual value of %g, primal value of %g\n",
     1305  // safeObjectiveValue,initialValue);
     1306  //safeObjectiveValue=max(2.0,1.0e-1*safeObjectiveValue);
     1307  //initialValue=max(100.0,1.0e-1*initialValue);;
     1308  //printf("temp safe dual value of %g, primal value of %g\n",
     1309  // safeObjectiveValue,initialValue);
    15291310  double zwLarge=1.0e2*initialValue;
    15301311  //zwLarge=1.0e40;
     
    15891370          if (objectiveValue>=0.0) {
    15901371            zVec_[iColumn]=objectiveValue + safeObjectiveValue*ratioZ;
    1591             wVec_[iColumn]=safeObjectiveValue*ratioW;
     1372            tVec_[iColumn]=safeObjectiveValue*ratioW;
    15921373          } else {
    15931374            zVec_[iColumn]=safeObjectiveValue*ratioZ;
    1594             wVec_[iColumn]=-objectiveValue + safeObjectiveValue*ratioW;
    1595           }
    1596           diagonal_[iColumn] = (t*s)/(s*wVec_[iColumn]+t*zVec_[iColumn]);
     1375            tVec_[iColumn]=-objectiveValue + safeObjectiveValue*ratioW;
     1376          }
     1377          diagonal_[iColumn] = (t*s)/(s*tVec_[iColumn]+t*zVec_[iColumn]);
    15971378        } else {
    15981379          //just lower bound
     
    16181399          if (objectiveValue>=0.0) {
    16191400            zVec_[iColumn]=objectiveValue + safeObjectiveValue*ratioZ;
    1620             wVec_[iColumn]=0.0;
     1401            tVec_[iColumn]=0.0;
    16211402          } else {
    16221403            zVec_[iColumn]=safeObjectiveValue*ratioZ;
    1623             wVec_[iColumn]=0.0;
     1404            tVec_[iColumn]=0.0;
    16241405          }
    16251406          diagonal_[iColumn] = s/zVec_[iColumn];
     
    16491430          if (objectiveValue>=0.0) {
    16501431            zVec_[iColumn]=0.0;
    1651             wVec_[iColumn]=safeObjectiveValue*ratioW;
     1432            tVec_[iColumn]=safeObjectiveValue*ratioW;
    16521433          } else {
    16531434            zVec_[iColumn]=0.0;
    1654             wVec_[iColumn]=-objectiveValue + safeObjectiveValue*ratioW;
    1655           }
    1656           diagonal_[iColumn] =  t/wVec_[iColumn];
     1435            tVec_[iColumn]=-objectiveValue + safeObjectiveValue*ratioW;
     1436          }
     1437          diagonal_[iColumn] =  t/tVec_[iColumn];
    16571438        } else {
    16581439          //free
    16591440          zVec_[iColumn]=safeObjectiveValue;
    1660           wVec_[iColumn]=safeObjectiveValue;
     1441          tVec_[iColumn]=safeObjectiveValue;
    16611442          newValue=fakeSolution[iColumn];
    16621443          if (newValue>=0.0) {
     
    16701451          }
    16711452          if (fabs(newValue)>1.0) {
    1672             diagonal_[iColumn]=fabs(newValue)/(zVec_[iColumn]+wVec_[iColumn]);
     1453            diagonal_[iColumn]=fabs(newValue)/(zVec_[iColumn]+tVec_[iColumn]);
    16731454          } else {
    1674             diagonal_[iColumn]=1.0/(zVec_[iColumn]+wVec_[iColumn]);
     1455            diagonal_[iColumn]=1.0/(zVec_[iColumn]+tVec_[iColumn]);
    16751456          }
    16761457          low=0.0;
     
    16861467      solution_[iColumn]=lower_[iColumn];
    16871468      zVec_[iColumn]=0.0;
    1688       wVec_[iColumn]=0.0;
     1469      tVec_[iColumn]=0.0;
    16891470      diagonal_[iColumn]=0.0;
    16901471    }
     
    16961477//phase 0=as is , 1 = after predictor , 2 after corrector
    16971478double ClpPredictorCorrector::complementarityGap(int & numberComplementarityPairs,
    1698                           const int phase)
     1479                                                 int & numberComplementarityItems,
     1480                                                 const int phase)
    16991481{
    17001482  double gap=0.0;
    17011483  //seems to be same coding for phase = 1 or 2
    17021484  numberComplementarityPairs=0;
     1485  numberComplementarityItems=0;
    17031486  double toleranceGap=0.0;
    17041487  double largestGap=0.0;
     
    17161499  dualTolerance=dualTolerance/scaleFactor_;
    17171500  double * zVec = zVec_;
    1718   double * wVec = wVec_;
     1501  double * tVec = tVec_;
    17191502  double * primal = solution_;
    17201503  double * lower = lower_;
     
    17221505  double * lowerSlack = lowerSlack_;
    17231506  double * upperSlack = upperSlack_;
    1724   double * work1 = deltaZ_;
    1725   double * work2 = deltaW_;
    1726   double * weights = weights_;
     1507  double * deltaZ = deltaZ_;
     1508  double * deltaT = deltaT_;
     1509  double * deltaX = deltaX_;
    17271510  for (int iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
    17281511    if (!fixedOrFree(iColumn)) {
    17291512      numberComplementarityPairs++;
    1730       //can collapse as if no lower bound both zVec and work1 0.0
     1513      //can collapse as if no lower bound both zVec and deltaZ 0.0
    17311514      if (lowerBound(iColumn)) {
     1515        numberComplementarityItems++;
    17321516        double dualValue;
    17331517        double primalValue;
     
    17381522          double change;
    17391523          if (!fakeLower(iColumn)) {
    1740             change =primal[iColumn]+weights[iColumn]-lowerSlack[iColumn]-lower[iColumn];
     1524            change =primal[iColumn]+deltaX[iColumn]-lowerSlack[iColumn]-lower[iColumn];
    17411525          } else {
    1742             change =weights[iColumn];
    1743           }
    1744           dualValue=zVec[iColumn]+actualDualStep_*work1[iColumn];
     1526            change =deltaX[iColumn];
     1527          }
     1528          dualValue=zVec[iColumn]+actualDualStep_*deltaZ[iColumn];
    17451529          primalValue=lowerSlack[iColumn]+actualPrimalStep_*change;
    17461530        }
     
    17671551      }
    17681552      if (upperBound(iColumn)) {
     1553        numberComplementarityItems++;
    17691554        double dualValue;
    17701555        double primalValue;
    17711556        if (!phase) {
    1772           dualValue=wVec[iColumn];
     1557          dualValue=tVec[iColumn];
    17731558          primalValue=upperSlack[iColumn];
    17741559        } else {
    17751560          double change;
    17761561          if (!fakeUpper(iColumn)) {
    1777             change =upper[iColumn]-primal[iColumn]-weights[iColumn]-upperSlack[iColumn];
     1562            change =upper[iColumn]-primal[iColumn]-deltaX[iColumn]-upperSlack[iColumn];
    17781563          } else {
    1779             change =-weights[iColumn];
    1780           }
    1781           dualValue=wVec[iColumn]+actualDualStep_*work2[iColumn];
     1564            change =-deltaX[iColumn];
     1565          }
     1566          dualValue=tVec[iColumn]+actualDualStep_*deltaT[iColumn];
    17821567          primalValue=upperSlack[iColumn]+actualPrimalStep_*change;
    17831568        }
     
    18241609  double extra =eExtra;
    18251610  double * zVec = zVec_;
    1826   double * wVec = wVec_;
     1611  double * tVec = tVec_;
    18271612  double * primal = solution_;
    1828   double * dual = dj_;
     1613  double * dj = dj_;
    18291614  double * lower = lower_;
    18301615  double * upper = upper_;
    18311616  double * lowerSlack = lowerSlack_;
    18321617  double * upperSlack = upperSlack_;
    1833   double * work2 = deltaW_;
    1834   double * work3 = deltaS_;
    1835   double * work4 = deltaT_;
     1618  double * workArray = workArray_;
    18361619
    18371620  int iColumn;
    1838   //printf("phase %d in setupForSolve, mu %g\n",phase,mu_);
     1621#ifdef SOME_DEBUG
     1622  printf("phase %d in setupForSolve, mu %g\n",phase,mu_);
     1623#endif
    18391624  switch (phase) {
    18401625  case 0:
     1626    CoinMemcpyN(errorRegion_,numberRows_,rhsB_);
    18411627    for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
     1628      rhsC_[iColumn]=0.0;
     1629      rhsU_[iColumn]=0.0;
     1630      rhsL_[iColumn]=0.0;
     1631      rhsZ_[iColumn]=0.0;
     1632      rhsT_[iColumn]=0.0;
    18421633      if (!flagged(iColumn)) {
    1843         double value= dual[iColumn];
     1634        rhsC_[iColumn] = dj[iColumn]-zVec[iColumn]+tVec[iColumn];
    18441635        if (lowerBound(iColumn)) {
    18451636          if (!fakeLower(iColumn)) {
    1846             value+=zVec[iColumn]*
    1847               (primal[iColumn]-lowerSlack[iColumn]-lower[iColumn])/
    1848                  (lowerSlack[iColumn]+extra);
     1637            rhsZ_[iColumn] = -zVec[iColumn]*(lowerSlack[iColumn]+extra);
     1638            rhsL_[iColumn] = primal[iColumn]-lowerSlack[iColumn]-lower[iColumn];
    18491639          }
    18501640        }
    18511641        if (upperBound(iColumn)) {
    18521642          if (!fakeUpper(iColumn)) {
    1853             value+=wVec[iColumn]*
    1854               (primal[iColumn]+upperSlack[iColumn]-upper[iColumn])/
    1855                  (upperSlack[iColumn]+extra);
    1856           }
    1857         }
    1858 #ifndef KKT
    1859         work2[iColumn]=diagonal_[iColumn]*value;
    1860 #else
    1861         work2[iColumn]=value;
    1862 #endif
    1863       } else {
    1864         work2[iColumn]=0.0;
    1865       }
     1643            rhsT_[iColumn] = -tVec[iColumn]*(upperSlack[iColumn]+extra);
     1644            rhsU_[iColumn] = -(primal[iColumn] + upperSlack[iColumn] - upper[iColumn]);
     1645          }
     1646        }
     1647      }
    18661648    }
    18671649    break;
    18681650  case 1:
     1651    // could be stored in delta2?
    18691652    for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
     1653      rhsC_[iColumn]=0.0;
     1654      //rhsU_[iColumn]=0.0;
     1655      //rhsL_[iColumn]=0.0;
     1656      rhsZ_[iColumn]=0.0;
     1657      rhsT_[iColumn]=0.0;
    18701658      if (!flagged(iColumn)) {
    1871         double value= 0.0;
     1659        rhsC_[iColumn] = dj[iColumn]-zVec[iColumn]+tVec[iColumn];
     1660        if (lowerBound(iColumn)) {
     1661          rhsZ_[iColumn] = mu_ -zVec[iColumn]*(lowerSlack[iColumn]+extra)
     1662            - deltaZ_[iColumn]*deltaX_[iColumn];
     1663          // To bring in line with OSL
     1664          //rhsZ_[iColumn] -= deltaZ_[iColumn]*rhsL_[iColumn];
     1665        }
     1666        if (upperBound(iColumn)) {
     1667          rhsT_[iColumn] = mu_ -tVec[iColumn]*(upperSlack[iColumn]+extra)
     1668            +deltaT_[iColumn]*deltaX_[iColumn];
     1669          // To bring in line with OSL
     1670          //rhsT_[iColumn] += deltaT_[iColumn]*rhsU_[iColumn];
     1671        }
     1672      }
     1673    }
     1674    break;
     1675  case 2:
     1676    CoinMemcpyN(errorRegion_,numberRows_,rhsB_);
     1677    for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
     1678      rhsC_[iColumn]=0.0;
     1679      rhsU_[iColumn]=0.0;
     1680      rhsL_[iColumn]=0.0;
     1681      rhsZ_[iColumn]=0.0;
     1682      rhsT_[iColumn]=0.0;
     1683      if (!flagged(iColumn)) {
     1684        rhsC_[iColumn] = dj[iColumn]-zVec[iColumn]+tVec[iColumn];
    18721685        if (lowerBound(iColumn)) {
    18731686          if (!fakeLower(iColumn)) {
    1874             value-=(mu_-work3[iColumn]-zVec[iColumn]*
    1875               (primal[iColumn]-lowerSlack[iColumn]-lower[iColumn]))/
    1876                  (lowerSlack[iColumn]+extra);
    1877           } else {
    1878             value-=(mu_-work3[iColumn])/(lowerSlack[iColumn]+extra);
     1687            rhsZ_[iColumn] = mu_ - zVec[iColumn]*(lowerSlack[iColumn]+extra);
     1688            rhsL_[iColumn] = primal[iColumn]-lowerSlack[iColumn]-lower[iColumn];
    18791689          }
    18801690        }
    18811691        if (upperBound(iColumn)) {
    18821692          if (!fakeUpper(iColumn)) {
    1883             value+=(mu_-work4[iColumn]+wVec[iColumn]*
    1884               (primal[iColumn]+upperSlack[iColumn]-upper[iColumn]))/
    1885                  (upperSlack[iColumn]+extra);
    1886           } else {
    1887             value+=(mu_-work4[iColumn])/(upperSlack[iColumn]+extra);
    1888           }
    1889         }
    1890 #ifndef KKT
    1891         work2[iColumn]=diagonal_[iColumn]*(dual[iColumn]+value);
    1892 #else
    1893         work2[iColumn]=dual[iColumn]+value;
    1894 #endif
    1895       } else {
    1896         work2[iColumn]=0.0;
    1897       }
    1898     }
    1899     break;
    1900   case 2:
    1901     for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
    1902       if (!flagged(iColumn)) {
    1903         double value= 0.0;
    1904         if (lowerBound(iColumn)) {
    1905           if (!fakeLower(iColumn)) {
    1906             value-=(mu_-zVec[iColumn]*
    1907               (primal[iColumn]-lowerSlack[iColumn]-lower[iColumn]))/
    1908                  (lowerSlack[iColumn]+extra);
    1909           } else {
    1910             value-=(mu_-zVec[iColumn])/ (lowerSlack[iColumn]+extra);
    1911           }
    1912         }
    1913         if (upperBound(iColumn)) {
    1914           if (!fakeUpper(iColumn)) {
    1915             value+=(mu_+wVec[iColumn]*
    1916               (primal[iColumn]+upperSlack[iColumn]-upper[iColumn]))/
    1917                  (upperSlack[iColumn]+extra);
    1918           } else {
    1919             value+=(mu_+wVec[iColumn])/ (upperSlack[iColumn]+extra);
    1920           }
    1921         }
    1922 #ifndef KKT
    1923         work2[iColumn]=diagonal_[iColumn]*(dual[iColumn]+value);
    1924 #else
    1925         work2[iColumn]=dual[iColumn]+value;
    1926 #endif
    1927       } else {
    1928         work2[iColumn]=0.0;
    1929       }
     1693            rhsT_[iColumn] = mu_ - tVec[iColumn]*(upperSlack[iColumn]+extra);
     1694            rhsU_[iColumn] = -(primal[iColumn] + upperSlack[iColumn] - upper[iColumn]);
     1695          }
     1696        }
     1697      }
    19301698    }
    19311699    break;
     
    19341702      double minBeta = 0.1*mu_;
    19351703      double maxBeta = 10.0*mu_;
    1936       double * work1 = deltaZ_;
    1937       double * weights = weights_;
     1704      double dualStep = min(1.0,actualDualStep_+0.1);
     1705      double primalStep = min(1.0,actualPrimalStep_+0.1);
     1706#ifdef SOME_DEBUG
     1707      printf("good complementarity range %g to %g\n",minBeta,maxBeta);
     1708#endif
     1709      //minBeta=0.0;
     1710      //maxBeta=COIN_DBL_MAX;
    19381711      for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
    1939         work3[iColumn]=0.0;
    1940         work4[iColumn]=0.0;
    19411712        if (!flagged(iColumn)) {
    1942           double value= 0.0;
    19431713          if (lowerBound(iColumn)) {
    1944             double change;
    1945             if (!fakeLower(iColumn)) {
    1946               change =primal[iColumn]+weights[iColumn]-lowerSlack[iColumn]-lower[iColumn];
    1947             } else {
    1948               change =weights[iColumn];
    1949             }
    1950             double dualValue=zVec[iColumn]+actualDualStep_*work1[iColumn];
    1951             double primalValue=lowerSlack[iColumn]+actualPrimalStep_*change;
     1714            double change = rhsL_[iColumn] + deltaX_[iColumn];
     1715            double dualValue=zVec[iColumn]+dualStep*deltaZ_[iColumn];
     1716            double primalValue=lowerSlack[iColumn]+primalStep*change;
    19521717            double gapProduct=dualValue*primalValue;
     1718            if (gapProduct>0.0&&dualValue<0.0)
     1719              gapProduct = - gapProduct;
     1720#ifdef FULL_DEBUG
     1721            delta2Z_[iColumn]=gapProduct;
     1722            if (delta2Z_[iColumn]<minBeta||delta2Z_[iColumn]>maxBeta)
     1723              printf("lower %d primal %g, dual %g, gap %g\n",
     1724                     iColumn,primalValue,dualValue,gapProduct);
     1725#endif
     1726            double value= 0.0;
    19531727            if (gapProduct<minBeta) {
    1954               value-= (minBeta-gapProduct)/(lowerSlack[iColumn]+extra);
     1728              value= 2.0*(minBeta-gapProduct);
     1729              value= (mu_-gapProduct);
     1730              value= (minBeta-gapProduct);
     1731              assert (value>0.0);
    19551732            } else if (gapProduct>maxBeta) {
    1956               value-= max(maxBeta-gapProduct,-maxBeta)/(lowerSlack[iColumn]+extra);
     1733              value= max(maxBeta-gapProduct,-maxBeta);
     1734              assert (value<0.0);
    19571735            }
    1958             value=-value;
    1959             work3[iColumn]=value;
     1736            //#define AGAIN
     1737#ifdef AGAIN
     1738            //rhsZ_[iColumn] = mu_ -zVec[iColumn]*(lowerSlack[iColumn]+extra)
     1739            //- deltaZ_[iColumn]*deltaX_[iColumn];
     1740            // To bring in line with OSL
     1741            rhsZ_[iColumn] += deltaZ_[iColumn]*rhsL_[iColumn];
     1742#endif
     1743            rhsZ_[iColumn] += value;
    19601744          } 
    19611745          if (upperBound(iColumn)) {
    1962             double change;
    1963             if (!fakeUpper(iColumn)) {
    1964               change =upper[iColumn]-primal[iColumn]-weights[iColumn]-upperSlack[iColumn];
    1965             } else {
    1966               change =-weights[iColumn];
    1967             }
    1968             double dualValue=wVec[iColumn]+actualDualStep_*work2[iColumn];
    1969             double primalValue=upperSlack[iColumn]+actualPrimalStep_*change;
     1746            double change = rhsU_[iColumn]-deltaX_[iColumn];
     1747            double dualValue=tVec[iColumn]+dualStep*deltaT_[iColumn];
     1748            double primalValue=upperSlack[iColumn]+primalStep*change;
    19701749            double gapProduct=dualValue*primalValue;
     1750            if (gapProduct>0.0&&dualValue<0.0)
     1751              gapProduct = - gapProduct;
     1752#ifdef FULL_DEBUG
     1753            delta2T_[iColumn]=gapProduct;
     1754            if (delta2T_[iColumn]<minBeta||delta2T_[iColumn]>maxBeta)
     1755              printf("upper %d primal %g, dual %g, gap %g\n",
     1756                     iColumn,primalValue,dualValue,gapProduct);
     1757#endif
     1758            double value= 0.0;
    19711759            if (gapProduct<minBeta) {
    1972               value+= (minBeta-gapProduct)/(upperSlack[iColumn]+extra);
     1760              value= (minBeta-gapProduct);
     1761              assert (value>0.0);
    19731762            } else if (gapProduct>maxBeta) {
    1974               value+= max(maxBeta-gapProduct,-maxBeta)/(upperSlack[iColumn]+extra);
     1763              value= max(maxBeta-gapProduct,-maxBeta);
     1764              assert (value<0.0);
    19751765            }
    1976             value=-value;
    1977             work4[iColumn]=value-work3[iColumn];
     1766#ifdef AGAIN
     1767            //rhsT_[iColumn] = mu_ -tVec[iColumn]*(upperSlack[iColumn]+extra)
     1768            //+deltaT_[iColumn]*deltaX_[iColumn];
     1769            // To bring in line with OSL
     1770            rhsT_[iColumn] -= deltaT_[iColumn]*rhsU_[iColumn];
     1771#endif
     1772            rhsT_[iColumn] += value;
    19781773          }
    1979 #ifndef KKT
    1980           work2[iColumn]=diagonal_[iColumn]*value;
    1981 #else
    1982           work2[iColumn]=value;
    1983 #endif
    1984         } else {
    1985           work2[iColumn]=0.0;
    19861774        }
    19871775      }
     
    19891777    break;
    19901778  } /* endswitch */
     1779#ifndef KKT
     1780  for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
     1781    double value = rhsC_[iColumn];
     1782    if (lowerBound(iColumn)) {
     1783      value -= rhsZ_[iColumn]/(lowerSlack[iColumn]+extra);
     1784      value += zVec[iColumn]*rhsL_[iColumn]/(lowerSlack[iColumn]+extra);
     1785    }
     1786    if (upperBound(iColumn)) {
     1787      value += rhsT_[iColumn]/(upperSlack[iColumn]+extra);
     1788      value -= tVec[iColumn]*rhsU_[iColumn]/(upperSlack[iColumn]+extra);
     1789    }
     1790    workArray[iColumn]=diagonal_[iColumn]*value;
     1791  }
     1792#else
     1793  for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
     1794    double value = rhsC_[iColumn];
     1795    if (lowerBound(iColumn)) {
     1796      value -= rhsZ_[iColumn]/(lowerSlack[iColumn]+extra);
     1797      value += zVec[iColumn]*rhsL_[iColumn]/(lowerSlack[iColumn]+extra);
     1798    }
     1799    if (upperBound(iColumn)) {
     1800      value += rhsT_[iColumn]/(upperSlack[iColumn]+extra);
     1801      value -= tVec[iColumn]*rhsU_[iColumn]/(upperSlack[iColumn]+extra);
     1802    }
     1803      workArray[iColumn]=value;
     1804  }
     1805#endif
    19911806}
    19921807//method: sees if looks plausible change in complementarity
     
    19961811  bool goodMove=false;
    19971812  int nextNumber;
     1813  int nextNumberItems;
    19981814  int numberTotal = numberRows_+numberColumns_;
    19991815  double returnGap=bestNextGap;
    2000   double nextGap=complementarityGap(nextNumber,2);
    2001   if (nextGap>bestNextGap)
     1816  double nextGap=complementarityGap(nextNumber,nextNumberItems,2);
     1817  if (nextGap>bestNextGap&&nextGap>-0.99*complementarityGap_) {
     1818#ifdef SOME_DEBUG
     1819    printf("checkGood phase 1 next gap %g, phase 0 %g, old gap %g\n",
     1820           nextGap,bestNextGap,complementarityGap_);
     1821#endif
    20021822    return false;
    2003   else
     1823  } else {
    20041824    returnGap=nextGap;
     1825  }
    20051826  double step;
    20061827  if (actualDualStep_>actualPrimalStep_) {
     
    20721893    double deltaObjectivePrimal=0.0;
    20731894    double deltaObjectiveDual=
    2074       innerProduct(updateRegion_,numberRows_,
     1895      innerProduct(deltaY_,numberRows_,
    20751896                   rhsFixRegion_);
    20761897    double error=0.0;
    2077     double * work3 = deltaS_;
    2078     CoinZeroN(work3,numberColumns_);
    2079     CoinMemcpyN(updateRegion_,numberRows_,work3+numberColumns_);
    2080     matrix_->transposeTimes(-1.0,updateRegion_,work3);
    2081     double * work1 = deltaZ_;
    2082     double * work2 = deltaW_;
     1898    double * workArray = workArray_;
     1899    CoinZeroN(workArray,numberColumns_);
     1900    CoinMemcpyN(deltaY_,numberRows_,workArray+numberColumns_);
     1901    matrix_->transposeTimes(-1.0,deltaY_,workArray);
     1902    double * deltaZ = deltaZ_;
     1903    double * deltaT = deltaT_;
    20831904    double * lower = lower_;
    20841905    double * upper = upper_;
    2085     //direction vector in weights
    2086     double * weights = weights_;
     1906    //direction vector in deltaX
     1907    double * deltaX = deltaX_;
    20871908    double * cost = cost_;
    20881909    //double sumPerturbCost=0.0;
     
    20901911      if (!flagged(iColumn)) {
    20911912        if (lowerBound(iColumn)) {
    2092           //sumPerturbCost+=weights[iColumn];
    2093           deltaObjectiveDual+=work1[iColumn]*lower[iColumn];
     1913          //sumPerturbCost+=deltaX[iColumn];
     1914          deltaObjectiveDual+=deltaZ[iColumn]*lower[iColumn];
    20941915        }
    20951916        if (upperBound(iColumn)) {
    2096           //sumPerturbCost-=weights[iColumn];
    2097           deltaObjectiveDual-=work2[iColumn]*upper[iColumn];
     1917          //sumPerturbCost-=deltaX[iColumn];
     1918          deltaObjectiveDual-=deltaT[iColumn]*upper[iColumn];
    20981919        }
    2099         double change = fabs(work3[iColumn]-work1[iColumn]+work2[iColumn]);
     1920        double change = fabs(workArray[iColumn]-deltaZ[iColumn]+deltaT[iColumn]);
    21001921        error = max (change,error);
    21011922      }
    2102       deltaObjectivePrimal += cost[iColumn] * weights[iColumn];
     1923      deltaObjectivePrimal += cost[iColumn] * deltaX[iColumn];
    21031924    }
    21041925    //deltaObjectivePrimal+=sumPerturbCost*linearPerturbation_;
     
    21471968  const double gammad = 1.0e-8;
    21481969  int nextNumber;
    2149   double nextGap=complementarityGap(nextNumber,2);
     1970  int nextNumberItems;
     1971  double nextGap=complementarityGap(nextNumber,nextNumberItems,2);
    21501972  if (nextGap>bestNextGap)
    21511973    return false;
     
    21531975  bool goodMove=true;
    21541976  double * deltaZ = deltaZ_;
    2155   double * deltaW = deltaW_;
    2156   double * deltaS = deltaS_;
    21571977  double * deltaT = deltaT_;
     1978  double * deltaSL = deltaSL_;
     1979  double * deltaSU = deltaSU_;
    21581980  double * zVec = zVec_;
    2159   double * wVec = wVec_;
     1981  double * tVec = tVec_;
    21601982  double * lowerSlack = lowerSlack_;
    21611983  double * upperSlack = upperSlack_;
     
    21631985    if (!flagged(iColumn)) {
    21641986      if (lowerBound(iColumn)) {
    2165         double part1=lowerSlack[iColumn]+actualPrimalStep_*deltaS[iColumn];
     1987        double part1=lowerSlack[iColumn]+actualPrimalStep_*deltaSL[iColumn];
    21661988        double part2=zVec[iColumn]+actualDualStep_*deltaZ[iColumn];
    21671989        if (part1*part2<lowerBoundGap) {
     
    21721994      if (upperBound(iColumn)) {
    21731995        double part1=upperSlack[iColumn]+actualPrimalStep_*deltaT[iColumn];
    2174         double part2=wVec[iColumn]+actualDualStep_*deltaW[iColumn];
     1996        double part2=tVec[iColumn]+actualDualStep_*deltaSU[iColumn];
    21751997        if (part1*part2<lowerBoundGap) {
    21761998          goodMove=false;
     
    22102032{
    22112033  //update pi
    2212   multiplyAdd(updateRegion_,numberRows_,actualDualStep_,dual_,1.0);
     2034  multiplyAdd(deltaY_,numberRows_,actualDualStep_,dual_,1.0);
    22132035  CoinZeroN(errorRegion_,numberRows_);
    22142036  CoinZeroN(rhsFixRegion_,numberRows_);
     
    22592081    qDiagonal=1.0e-8*mu_;
    22602082  }
    2261   double norm=1.0e-12;
    22622083  double widenGap=1.0e1;
    22632084  //largest allowable ratio of lowerSlack/zVec (etc)
     
    22812102  double smallestDiagonal=1.0e50;
    22822103  double * zVec = zVec_;
    2283   double * wVec = wVec_;
     2104  double * tVec = tVec_;
    22842105  double * primal = solution_;
    22852106  double * dual = dj_;
     
    22882109  double * lowerSlack = lowerSlack_;
    22892110  double * upperSlack = upperSlack_;
    2290   double * work1 = deltaZ_;
    2291   double * work2 = deltaW_;
     2111  double * deltaZ = deltaZ_;
     2112  double * deltaT = deltaT_;
    22922113  double * diagonal = diagonal_;
    2293   //direction vector in weights
    2294   double * weights = weights_;
     2114  //direction vector in deltaX
     2115  double * deltaX = deltaX_;
    22952116  double * cost = cost_;
    22962117  for (int iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
     
    22982119      double reducedCost=dual[iColumn];
    22992120      bool thisKilled=false;
    2300       double zValue = zVec[iColumn] + actualDualStep_*work1[iColumn];
    2301       double wValue = wVec[iColumn] + actualDualStep_*work2[iColumn];
     2121      double zValue = zVec[iColumn] + actualDualStep_*deltaZ[iColumn];
     2122      double tValue = tVec[iColumn] + actualDualStep_*deltaT[iColumn];
    23022123      zVec[iColumn]=zValue;
    2303       wVec[iColumn]=wValue;
    2304       double thisWeight=weights[iColumn];
     2124      tVec[iColumn]=tValue;
     2125      double thisWeight=deltaX[iColumn];
    23052126      double oldPrimal = primal[iColumn];
    23062127      double newPrimal=primal[iColumn]+actualPrimalStep_*thisWeight;
     
    23212142      double smallerSlack=COIN_DBL_MAX;
    23222143      double largerzw;
    2323       if (zValue>wValue) {
     2144      if (zValue>tValue) {
    23242145        largerzw=zValue;
    23252146      } else {
    2326         largerzw=wValue;
     2147        largerzw=tValue;
    23272148      }
    23282149      if (lowerBound(iColumn)) {
     
    24412262        }
    24422263        //for changing slack
    2443         double wValue2 = wValue;
    2444         if (wValue2<epsilonBase) {
    2445           wValue2=epsilonBase;
     2264        double tValue2 = tValue;
     2265        if (tValue2<epsilonBase) {
     2266          tValue2=epsilonBase;
    24462267        }
    24472268        //make sure reasonable
    2448         if (wValue<epsilon) {
    2449           wValue=epsilon;
    2450         }
    2451         //store modified wVec
    2452         //no wVec[iColumn]=wValue;
     2269        if (tValue<epsilon) {
     2270          tValue=epsilon;
     2271        }
     2272        //store modified tVec
     2273        //no tVec[iColumn]=tValue;
    24532274        double feasibleSlack=upper[iColumn]-newPrimal;
    24542275        if (feasibleSlack>0.0&&newSlack>0.0) {
     
    24682289            }
    24692290            //set FAKE here
    2470             if (newSlack>wValue2*largestRatio&&newSlack>smallGap2) {
     2291            if (newSlack>tValue2*largestRatio&&newSlack>smallGap2) {
    24712292              setFakeUpper(iColumn);
    2472               newSlack=wValue2*largestRatio;
     2293              newSlack=tValue2*largestRatio;
    24732294              if (newSlack<smallGap2) {
    24742295                newSlack=smallGap2;
     
    24772298            }
    24782299          } else {
    2479             newSlack=wValue2*largestRatio;
     2300            newSlack=tValue2*largestRatio;
    24802301            if (newSlack<smallGap2) {
    24812302              newSlack=smallGap2;
     
    24952316          }
    24962317        }
    2497         if (wVec[iColumn]>dualTolerance) {
    2498           dualObjectiveThis-=upper[iColumn]*wVec[iColumn];
     2318        if (tVec[iColumn]>dualTolerance) {
     2319          dualObjectiveThis-=upper[iColumn]*tVec[iColumn];
    24992320        }
    25002321        upperSlack[iColumn]=newSlack;
     
    25262347        }
    25272348        zVec[iColumn]=gap*freeMultiplier;
    2528         wVec[iColumn]=zVec[iColumn];
     2349        tVec[iColumn]=zVec[iColumn];
    25292350        //fake to give correct result
    25302351        s=1.0;
    25312352        t=1.0;
    2532         wValue=0.0;
     2353        tValue=0.0;
    25332354        zValue=2.0*freeMultiplier+qDiagonal;
    25342355      }
    25352356      primal[iColumn]=newPrimal;
    2536       if (fabs(newPrimal)>norm) {
    2537         norm=fabs(newPrimal);
     2357      if (fabs(newPrimal)>solutionNorm) {
     2358        solutionNorm=fabs(newPrimal);
    25382359      }
    25392360      if (!thisKilled) {
     
    25462367          t=largeGap;
    25472368        }
    2548         double divisor = s*wValue+t*zValue;
     2369        double divisor = s*tValue+t*zValue;
    25492370        double diagonalValue=(t*s)/divisor;
    25502371        diagonal[iColumn]=diagonalValue;
     
    25632384          smallestDiagonal=diagonalValue;
    25642385        }
    2565         double dualInfeasibility=reducedCost-zVec[iColumn]+wVec[iColumn];
     2386        double dualInfeasibility=reducedCost-zVec[iColumn]+tVec[iColumn];
    25662387        if (fabs(dualInfeasibility)>dualTolerance) {
    25672388          if ((!lowerBound(iColumn))&&
     
    25762397          maximumDualError=dualInfeasibility;
    25772398        }
    2578         work1[iColumn]=0.0;
     2399        deltaZ[iColumn]=0.0;
    25792400      } else {
    25802401        numberKilled++;
    25812402        diagonal[iColumn]=0.0;
    25822403        zVec[iColumn]=0.0;
    2583         wVec[iColumn]=0.0;
     2404        tVec[iColumn]=0.0;
    25842405        setFlagged(iColumn);
    25852406        setFixedOrFree(iColumn);
    2586         work1[iColumn]=newPrimal;
     2407        deltaZ[iColumn]=newPrimal;
    25872408        offsetObjective+=newPrimal*cost[iColumn];
    25882409      }
    25892410    } else {
    2590       work1[iColumn]=primal[iColumn];
     2411      deltaZ[iColumn]=primal[iColumn];
    25912412      diagonal[iColumn]=0.0;
    25922413      offsetObjective+=primal[iColumn]*cost[iColumn];
     
    26192440#endif
    26202441  // update rhs region
    2621   multiplyAdd(work1+numberColumns_,numberRows_,-1.0,rhsFixRegion_,1.0);
    2622   matrix_->times(1.0,work1,rhsFixRegion_);
     2442  multiplyAdd(deltaZ+numberColumns_,numberRows_,-1.0,rhsFixRegion_,1.0);
     2443  matrix_->times(1.0,deltaZ,rhsFixRegion_);
    26232444  primalObjectiveValue=primalObjectiveValue+offsetObjective;
    26242445  dualObjectiveValue+=offsetObjective+dualFake;
     
    26372458  sumDualInfeasibilities_=maximumDualError;
    26382459  maximumBoundInfeasibility_ = maximumBoundInfeasibility;
    2639   solutionNorm_ = norm;
    26402460  //compute error and fixed RHS
    26412461  multiplyAdd(solution_+numberColumns_,numberRows_,-1.0,errorRegion_,0.0);
     
    27542574ClpPredictorCorrector::affineProduct()
    27552575{
    2756   double * work1 = deltaZ_;
    2757   double * work2 = deltaW_;
    2758   double * work3 = deltaS_;
    2759   double * work4 = deltaT_;
     2576  double * deltaZ = deltaZ_;
     2577  double * deltaT = deltaT_;
    27602578  double * lower = lower_;
    27612579  double * upper = upper_;
     
    27632581  double * upperSlack = upperSlack_;
    27642582  double * primal = solution_;
    2765   //direction vector in weights
    2766   double * weights = weights_;
     2583  //direction vector in deltaX
     2584  double * deltaX = deltaX_;
    27672585  double product = 0.0;
    2768   //IF zVec starts as 0 then work1 always zero
     2586  //IF zVec starts as 0 then deltaZ always zero
    27692587  //(remember if free then zVec not 0)
    27702588  //I think free can be done with careful use of boundSlacks to zero
    27712589  //out all we want
    27722590  for (int iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
    2773     double w3=work1[iColumn]*weights[iColumn];
    2774     double w4=-work2[iColumn]*weights[iColumn];
     2591    double w3=deltaZ[iColumn]*deltaX[iColumn];
     2592    double w4=-deltaT[iColumn]*deltaX[iColumn];
    27752593    if (lowerBound(iColumn)) {
    27762594      if (!fakeLower(iColumn)) {
    2777         w3+=work1[iColumn]*(primal[iColumn]-lowerSlack[iColumn]-lower[iColumn]);
     2595        w3+=deltaZ[iColumn]*(primal[iColumn]-lowerSlack[iColumn]-lower[iColumn]);
    27782596      }
    27792597      product+=w3;
     
    27812599    if (upperBound(iColumn)) {
    27822600      if (!fakeUpper(iColumn)) {
    2783         w4+=work2[iColumn]*(-primal[iColumn]-upperSlack[iColumn]+upper[iColumn]);
     2601        w4+=deltaT[iColumn]*(-primal[iColumn]-upperSlack[iColumn]+upper[iColumn]);
    27842602      }
    27852603      product+=w4;
    27862604    }
    2787     work3[iColumn]=w3;
    2788     work4[iColumn]=w4;
    27892605  }
    27902606  return product;
  • trunk/ClpSolve.cpp

    r340 r369  
    892892    model2->setNumberIterations(model2->numberIterations()+totalIterations);
    893893  } else if (method==ClpSolve::useBarrier) {
    894     printf("***** experimental pretty crude barrier\n");
     894    //printf("***** experimental pretty crude barrier\n");
    895895    ClpInterior barrier(*model2);
    896896    if (interrupt)
     
    898898#ifdef REAL_BARRIER
    899899    // uncomment this if you have Anshul Gupta's wsmp package
    900     //ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp(max(100,model2->numberRows()/10));
    901     ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp();
     900    ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp(max(100,model2->numberRows()/10));
     901    //ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp();
    902902    barrier.setCholesky(cholesky);
    903903    //ClpCholeskyWssmpKKT * cholesky = new ClpCholeskyWssmpKKT(max(100,model2->numberRows()/10));
     
    910910#ifndef SAVEIT
    911911    barrier.primalDual();
     912    //printf("********** Stopping as this is debug run\n");
     913    //exit(99);
    912914#elif SAVEIT==1
    913915    barrier.primalDual();
     
    930932      <<CoinMessageEol;
    931933    timeX=time2;
    932     if (barrier.maximumBarrierIterations()&&barrier.status()!=4) {
     934    if (barrier.maximumBarrierIterations()&&barrier.status()<4) {
    933935      printf("***** crossover - needs more thought on difficult models\n");
    934936      // move solutions
  • trunk/include/ClpInterior.hpp

    r268 r369  
    77   John Tomlin (pdco)
    88   John Forrest (standard predictor-corrector)
     9
     10   Note JJF has added arrays - this takes more memory but makes
     11   flow easier to understand and hopefully easier to extend
    912
    1013 */
     
    418421  /// rhsFixRegion.
    419422  double * rhsFixRegion_;
    420   /// updateRegion.
    421   double * updateRegion_;
    422423  /// upperSlack
    423424  double * upperSlack_;
     
    426427  /// diagonal
    427428  double * diagonal_;
    428   /// weights
    429   double * weights_;
    430429  /// solution
    431430  double * solution_;
    432   /// work1 or deltaZ.
     431  /// work array
     432  double * workArray_;
     433  /// delta X
     434  double * deltaX_;
     435  /// delta Y
     436  double * deltaY_;
     437  /// deltaZ.
    433438  double * deltaZ_;
    434   /// work2 or deltaW.
    435   double * deltaW_;
    436   /// work3 or deltaS.
    437   double * deltaS_;
    438   /// work4 or deltaT.
     439  /// deltaT.
    439440  double * deltaT_;
     441  /// deltaS.
     442  double * deltaSU_;
     443  double * deltaSL_;
     444  /// rhs B
     445  double * rhsB_;
     446  /// rhsU.
     447  double * rhsU_;
     448  /// rhsL.
     449  double * rhsL_;
     450  /// rhsZ.
     451  double * rhsZ_;
     452  /// rhsT.
     453  double * rhsT_;
     454  /// rhs C
     455  double * rhsC_;
    440456  /// zVec
    441457  double * zVec_;
    442   /// wVec
    443   double * wVec_;
     458  /// tVec
     459  double * tVec_;
    444460  /// cholesky.
    445461  ClpCholeskyBase * cholesky_;
    446   /// numberComplementarityPairs;
     462  /// numberComplementarityPairs i.e. ones with lower and/or upper bounds (not fixed)
    447463  int numberComplementarityPairs_;
     464  /// numberComplementarityItems_ i.e. number of active bounds
     465  int numberComplementarityItems_;
    448466  /// Maximum iterations
    449467  int maximumBarrierIterations_;
  • trunk/include/ClpPredictorCorrector.hpp

    r364 r369  
    4444  //         1 corrector
    4545  //         2 primal dual
    46   double findStepLength(const int phase,const double * oldWeight);
     46  double findStepLength( int phase);
    4747  /// findDirectionVector.
    4848  double findDirectionVector(const int phase);
     
    5151  /// complementarityGap.  Computes gap
    5252  //phase 0=as is , 1 = after predictor , 2 after corrector
    53   double complementarityGap(int & numberComplementarityPairs,
     53  double complementarityGap(int & numberComplementarityPairs,int & numberComplementarityItems,
    5454                            const int phase);
    5555  /// setupForSolve.
Note: See TracChangeset for help on using the changeset viewer.