source: trunk/Clp/src/ClpPredictorCorrector.cpp @ 1366

Last change on this file since 1366 was 1366, checked in by forrest, 11 years ago

try and improve quadratic and barrier

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 125.3 KB
Line 
1// Copyright (C) 2003, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4/*
5   Implements crude primal dual predictor corrector algorithm
6
7 */
8//#define SOME_DEBUG
9
10#include "CoinPragma.hpp"
11#include <math.h>
12
13#include "CoinHelperFunctions.hpp"
14#include "ClpPredictorCorrector.hpp"
15#include "CoinPackedMatrix.hpp"
16#include "ClpMessage.hpp"
17#include "ClpCholeskyBase.hpp"
18#include "ClpHelperFunctions.hpp"
19#include "ClpQuadraticObjective.hpp"
20#include <cfloat>
21#include <cassert>
22#include <string>
23#include <cstdio>
24#include <iostream>
25#if 0
26static int yyyyyy=0;
27void ClpPredictorCorrector::saveSolution(std::string fileName)
28{
29  FILE * fp=fopen(fileName.c_str(),"wb");
30  if (fp) {
31    int numberRows=numberRows_;
32    int numberColumns=numberColumns_;
33    fwrite(&numberRows,sizeof(int),1,fp);
34    fwrite(&numberColumns,sizeof(int),1,fp);
35    double dsave[20];
36    memset(dsave,0,sizeof(dsave));
37    fwrite(dsave,sizeof(double),20,fp);
38    int msave[20];
39    memset(msave,0,sizeof(msave));
40    msave[0]=numberIterations_;
41    fwrite(msave,sizeof(int),20,fp);
42    fwrite(dual_,sizeof(double),numberRows,fp);
43    fwrite(errorRegion_,sizeof(double),numberRows,fp);
44    fwrite(rhsFixRegion_,sizeof(double),numberRows,fp);
45    fwrite(solution_,sizeof(double),numberColumns,fp);
46    fwrite(solution_+numberColumns,sizeof(double),numberRows,fp);
47    fwrite(diagonal_,sizeof(double),numberColumns,fp);
48    fwrite(diagonal_+numberColumns,sizeof(double),numberRows,fp);
49    fwrite(wVec_,sizeof(double),numberColumns,fp);
50    fwrite(wVec_+numberColumns,sizeof(double),numberRows,fp);
51    fwrite(zVec_,sizeof(double),numberColumns,fp);
52    fwrite(zVec_+numberColumns,sizeof(double),numberRows,fp);
53    fwrite(upperSlack_,sizeof(double),numberColumns,fp);
54    fwrite(upperSlack_+numberColumns,sizeof(double),numberRows,fp);
55    fwrite(lowerSlack_,sizeof(double),numberColumns,fp);
56    fwrite(lowerSlack_+numberColumns,sizeof(double),numberRows,fp);
57    fclose(fp);
58  } else {
59    std::cout<<"Unable to open file "<<fileName<<std::endl;
60  }
61}
62#endif
63static double eScale=1.0e27;
64static double eBaseCaution=1.0e-12;
65static double eBase=1.0e-12;
66static double eRatio=1.0e40;
67static double eRatioCaution=1.0e25;
68static double eDiagonal=1.0e25;
69static double eDiagonalCaution=1.0e18;
70static double eExtra=1.0e-12;
71
72// main function
73
74int ClpPredictorCorrector::solve ( )
75{
76  problemStatus_=-1;
77  algorithm_=1;
78  //create all regions
79  if (!createWorkingData()) {
80    problemStatus_=4;
81    return 2;
82  }
83  //diagonalPerturbation_=1.0e-25;
84  ClpMatrixBase * saveMatrix = NULL;
85  // If quadratic then make copy so we can actually scale or normalize
86#ifndef NO_RTTI
87  ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
88#else
89  ClpQuadraticObjective * quadraticObj = NULL;
90  if (objective_->type()==2)
91    quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
92#endif
93  /* If modeSwitch is :
94     0 - normal
95     1 - bit switch off centering
96     2 - bit always do type 2
97     4 - accept corrector nearly always
98  */
99  int modeSwitch=0;
100  //if (quadraticObj)
101    //modeSwitch |= 1; // switch off centring for now
102  //if (quadraticObj)
103  //modeSwitch |=4;
104  ClpObjective * saveObjective = NULL;
105  if (quadraticObj) {
106    // check KKT is on
107    if (!cholesky_->kkt()) {
108      //No!
109      handler_->message(CLP_BARRIER_KKT,messages_)
110        <<CoinMessageEol;
111      return -1;
112    }
113    saveObjective=objective_;
114    // We are going to make matrix full rather half
115    objective_ = new ClpQuadraticObjective(*quadraticObj,1);
116  }
117  bool allowIncreasingGap = (modeSwitch&4)!=0;
118  // If scaled then really scale matrix
119  if (scalingFlag_>0&&rowScale_) {
120    saveMatrix = matrix_;
121    matrix_ = matrix_->scaledColumnCopy(this);
122  }
123  //initializeFeasible(); - this just set fixed flag
124  smallestInfeasibility_=COIN_DBL_MAX;
125  int i;
126  for (i=0;i<LENGTH_HISTORY;i++) 
127    historyInfeasibility_[i]=COIN_DBL_MAX;
128
129  //bool firstTime=true;
130  //firstFactorization(true);
131  int returnCode = cholesky_->order(this);
132  if (returnCode||cholesky_->symbolic()) {
133    printf("Error return from symbolic - probably not enough memory\n");
134    problemStatus_=4;
135    //delete all temporary regions
136    deleteWorkingData();
137    if (saveMatrix) {
138      // restore normal copy
139      delete matrix_;
140      matrix_ = saveMatrix;
141    }
142    // Restore quadratic objective if necessary
143    if (saveObjective) {
144      delete objective_;
145      objective_=saveObjective;
146    }
147    return -1;
148  }
149  mu_=1.0e10;
150  diagonalScaleFactor_=1.0;
151  //set iterations
152  numberIterations_=-1;
153  int numberTotal = numberRows_+numberColumns_;
154  //initialize solution here
155  if(createSolution()<0) {
156    printf("Not enough memory\n");
157    problemStatus_=4;
158    //delete all temporary regions
159    deleteWorkingData();
160    if (saveMatrix) {
161      // restore normal copy
162      delete matrix_;
163      matrix_ = saveMatrix;
164    }
165    return -1;
166  }
167  // Could try centering steps without any original step i.e. just center
168  //firstFactorization(false);
169  CoinZeroN(dual_,numberRows_);
170  multiplyAdd(solution_+numberColumns_,numberRows_,-1.0,errorRegion_,0.0);
171  matrix_->times(1.0,solution_,errorRegion_);
172  maximumRHSError_=maximumAbsElement(errorRegion_,numberRows_);
173  maximumBoundInfeasibility_=maximumRHSError_;
174  //double maximumDualError_=COIN_DBL_MAX;
175  //initialize
176  actualDualStep_=0.0;
177  actualPrimalStep_=0.0;
178  gonePrimalFeasible_=false;
179  goneDualFeasible_=false;
180  //bool hadGoodSolution=false;
181  diagonalNorm_=solutionNorm_;
182  mu_=solutionNorm_;
183  int numberFixed=updateSolution(-COIN_DBL_MAX);
184  int numberFixedTotal=numberFixed;
185  //int numberRows_DroppedBefore=0;
186  //double extra=eExtra;
187  //double maximumPerturbation=COIN_DBL_MAX;
188  //constants for infeas interior point
189  const double beta2 = 0.99995;
190  const double tau   = 0.00002;
191  double lastComplementarityGap=COIN_DBL_MAX;
192  double lastStep=1.0;
193  // use to see if to take affine
194  double checkGap = COIN_DBL_MAX;
195  int lastGoodIteration=0;
196  double bestObjectiveGap=COIN_DBL_MAX;
197  double bestObjective=COIN_DBL_MAX;
198  int saveIteration=-1;
199  bool sloppyOptimal=false;
200  double * savePi=NULL;
201  double * savePrimal=NULL;
202  // Extra regions for centering
203  double * saveX = new double[numberTotal];
204  double * saveY = new double[numberRows_];
205  double * saveZ = new double[numberTotal];
206  double * saveW = new double[numberTotal];
207  double * saveSL = new double[numberTotal];
208  double * saveSU = new double[numberTotal];
209  // Save smallest mu used in primal dual moves
210  double smallestPrimalDualMu=COIN_DBL_MAX;
211  double objScale = optimizationDirection_/
212    (rhsScale_*objectiveScale_);
213  while (problemStatus_<0) {
214    //#define FULL_DEBUG
215#ifdef FULL_DEBUG
216    {
217      int i;
218      printf("row    pi          artvec       rhsfx\n");
219      for (i=0;i<numberRows_;i++) {
220        printf("%d %g %g %g\n",i,dual_[i],errorRegion_[i],rhsFixRegion_[i]);
221      }
222      printf(" col  dsol  ddiag  dwvec  dzvec dbdslu dbdsll\n");
223      for (i=0;i<numberColumns_+numberRows_;i++) {
224        printf(" %d %g %g %g %g %g %g\n",i,solution_[i],diagonal_[i],wVec_[i],
225               zVec_[i],upperSlack_[i],lowerSlack_[i]);
226      }
227    }
228#endif
229    complementarityGap_=complementarityGap(numberComplementarityPairs_,
230                                           numberComplementarityItems_,0);
231    handler_->message(CLP_BARRIER_ITERATION,messages_)
232      <<numberIterations_
233      <<primalObjective_*objScale- dblParam_[ClpObjOffset]
234      << dualObjective_*objScale- dblParam_[ClpObjOffset]
235      <<complementarityGap_
236      <<numberFixedTotal
237      <<cholesky_->rank()
238      <<CoinMessageEol;
239#if 0
240    if (numberIterations_==-1) {
241      saveSolution("xxx.sav");
242      if (yyyyyy)
243        exit(99);
244    }
245#endif
246    // move up history
247    for (i=1;i<LENGTH_HISTORY;i++) 
248      historyInfeasibility_[i-1]=historyInfeasibility_[i];
249    historyInfeasibility_[LENGTH_HISTORY-1]=complementarityGap_;
250    // switch off saved if changes
251    //if (saveIteration+10<numberIterations_&&
252    //complementarityGap_*2.0<historyInfeasibility_[0])
253    //saveIteration=-1;
254    lastStep = CoinMin(actualPrimalStep_,actualDualStep_);
255    double goodGapChange;
256    //#define KEEP_GOING_IF_FIXED 0
257#ifndef KEEP_GOING_IF_FIXED
258#define KEEP_GOING_IF_FIXED 10000
259#endif
260    if (!sloppyOptimal) {
261      goodGapChange=0.93;
262    } else {
263      goodGapChange=0.7;
264      if (numberFixed>KEEP_GOING_IF_FIXED)
265        goodGapChange=0.99; // make more likely to carry on
266    } 
267    double gapO;
268    double lastGood=bestObjectiveGap;
269    if (gonePrimalFeasible_&&goneDualFeasible_) {
270      double largestObjective;
271      if (fabs(primalObjective_)>fabs(dualObjective_)) {
272        largestObjective = fabs(primalObjective_);
273      } else {
274        largestObjective = fabs(dualObjective_);
275      } 
276      if (largestObjective<1.0) {
277        largestObjective=1.0;
278      } 
279      gapO=fabs(primalObjective_-dualObjective_)/largestObjective;
280      handler_->message(CLP_BARRIER_OBJECTIVE_GAP,messages_)
281        <<gapO
282        <<CoinMessageEol;
283      //start saving best
284      if (gapO<bestObjectiveGap) 
285        bestObjectiveGap=gapO;
286      if (primalObjective_<bestObjective) {
287        saveIteration=numberIterations_;
288        bestObjective=primalObjective_;
289        if (!savePi) {
290          savePi=new double[numberRows_];
291          savePrimal = new double [numberTotal];
292        } 
293        CoinMemcpyN(dual_,numberRows_,savePi);
294        CoinMemcpyN(solution_,numberTotal,savePrimal);
295      } else if(gapO>2.0*bestObjectiveGap) {
296        //maybe be more sophisticated e.g. re-initialize having
297        //fixed variables and dropped rows
298        //std::cout <<" gap increasing "<<std::endl;
299      } 
300      //std::cout <<"could stop"<<std::endl;
301      //gapO=0.0;
302      if (fabs(primalObjective_-dualObjective_)<dualTolerance()) {
303        gapO=0.0;
304      } 
305    } else {
306      gapO=COIN_DBL_MAX;
307      if (saveIteration>=0) {
308        handler_->message(CLP_BARRIER_GONE_INFEASIBLE,messages_)
309          <<CoinMessageEol;
310        if (sloppyOptimal) {
311          // vaguely optimal
312          double scaledRHSError=maximumRHSError_/solutionNorm_;
313          if (maximumBoundInfeasibility_>1.0e-2||
314              scaledRHSError>1.0e-2||
315              maximumDualError_>objectiveNorm_*1.0e-2) {
316            handler_->message(CLP_BARRIER_EXIT2,messages_)
317              <<saveIteration
318              <<CoinMessageEol;
319            problemStatus_=0; // benefit of doubt
320            break;
321          }
322        } else {
323          // not close to optimal but check if getting bad
324          double scaledRHSError=maximumRHSError_/solutionNorm_;
325          if ((maximumBoundInfeasibility_>1.0e-1||
326              scaledRHSError>1.0e-1||
327               maximumDualError_>objectiveNorm_*1.0e-1)
328              &&(numberIterations_>50
329                 &&complementarityGap_>0.9*historyInfeasibility_[0])) {
330            handler_->message(CLP_BARRIER_EXIT2,messages_)
331              <<saveIteration
332              <<CoinMessageEol;
333            break;
334          }
335          if (complementarityGap_>0.95*checkGap&&bestObjectiveGap<1.0e-3&&
336              (numberIterations_>saveIteration+5||numberIterations_>100)) {
337            handler_->message(CLP_BARRIER_EXIT2,messages_)
338              <<saveIteration
339              <<CoinMessageEol;
340            break;
341          }
342        }
343      } 
344      if (complementarityGap_>0.5*checkGap&&primalObjective_>
345          bestObjective+1.0e-9&&
346          (numberIterations_>saveIteration+5||numberIterations_>100)) {
347        handler_->message(CLP_BARRIER_EXIT2,messages_)
348          <<saveIteration
349          <<CoinMessageEol;
350        break;
351      }
352    } 
353    if ((gapO<1.0e-6||(gapO<1.0e-4&&complementarityGap_<0.1))&&!sloppyOptimal) {
354      sloppyOptimal=true;
355      handler_->message(CLP_BARRIER_CLOSE_TO_OPTIMAL,messages_)
356        <<numberIterations_<<complementarityGap_
357        <<CoinMessageEol;
358    } 
359    int numberBack = quadraticObj ? 10 : 5;
360    //tryJustPredictor=true;
361    //printf("trying just predictor\n");
362    //}
363    if (complementarityGap_>=1.05*lastComplementarityGap) {
364      handler_->message(CLP_BARRIER_COMPLEMENTARITY,messages_)
365        <<complementarityGap_<<"increasing"
366        <<CoinMessageEol;
367      if (saveIteration>=0&&sloppyOptimal) {
368        handler_->message(CLP_BARRIER_EXIT2,messages_)
369          <<saveIteration
370          <<CoinMessageEol;
371        break;
372      } else if (numberIterations_-lastGoodIteration>=numberBack&&
373                 complementarityGap_<1.0e-6) {
374        break; // not doing very well - give up
375      } 
376    } else if (complementarityGap_<goodGapChange*lastComplementarityGap) {
377      lastGoodIteration=numberIterations_;
378      lastComplementarityGap=complementarityGap_;
379    } else if (numberIterations_-lastGoodIteration>=numberBack&&
380               complementarityGap_<1.0e-3) {
381      handler_->message(CLP_BARRIER_COMPLEMENTARITY,messages_)
382        <<complementarityGap_<<"not decreasing"
383        <<CoinMessageEol;
384      if (gapO>0.75*lastGood&&numberFixed<KEEP_GOING_IF_FIXED) {
385        break;
386      } 
387    } else if (numberIterations_-lastGoodIteration>=2&&
388               complementarityGap_<1.0e-6) {
389      handler_->message(CLP_BARRIER_COMPLEMENTARITY,messages_)
390        <<complementarityGap_<<"not decreasing"
391        <<CoinMessageEol;
392      break;
393    } 
394    if (numberIterations_>maximumBarrierIterations_||hitMaximumIterations()) {
395      handler_->message(CLP_BARRIER_STOPPING,messages_)
396        <<CoinMessageEol;
397      problemStatus_=3;
398      onStopped(); // set secondary status
399      break;
400    } 
401    if (gapO<targetGap_) {
402      problemStatus_=0;
403      handler_->message(CLP_BARRIER_EXIT,messages_)
404        <<" "
405        <<CoinMessageEol;
406        break;//finished
407    } 
408    if (complementarityGap_<1.0e-12) {
409      problemStatus_=0;
410      handler_->message(CLP_BARRIER_EXIT,messages_)
411        <<"- small complementarity gap"
412        <<CoinMessageEol;
413        break;//finished
414    } 
415    if (complementarityGap_<1.0e-10&&gapO<1.0e-10) {
416      problemStatus_=0;
417      handler_->message(CLP_BARRIER_EXIT,messages_)
418        <<"- objective gap and complementarity gap both small"
419        <<CoinMessageEol;
420        break;//finished
421    } 
422    if (gapO<1.0e-9) {
423      double value=gapO*complementarityGap_;
424      value*=actualPrimalStep_;
425      value*=actualDualStep_;
426      //std::cout<<value<<std::endl;
427      if (value<1.0e-17&&numberIterations_>lastGoodIteration) {
428        problemStatus_=0;
429        handler_->message(CLP_BARRIER_EXIT,messages_)
430          <<"- objective gap and complementarity gap both smallish and small steps"
431          <<CoinMessageEol;
432        break;//finished
433      } 
434    } 
435    double nextGap=COIN_DBL_MAX;
436    int nextNumber=0;
437    int nextNumberItems=0;
438    worstDirectionAccuracy_=0.0;
439    int newDropped=0;
440    //Predictor step
441    //prepare for cholesky.  Set up scaled diagonal in deltaX
442    //  ** for efficiency may be better if scale factor known before
443    double norm2=0.0;
444    double maximumValue;
445    getNorms(diagonal_,numberTotal,maximumValue,norm2);
446    diagonalNorm_ = sqrt(norm2/numberComplementarityPairs_);
447    diagonalScaleFactor_=1.0;
448    double maximumAllowable=eScale;
449    //scale so largest is less than allowable ? could do better
450    double factor=0.5;
451    while (maximumValue>maximumAllowable) {
452      diagonalScaleFactor_*=factor;
453      maximumValue*=factor;
454    } /* endwhile */
455    if (diagonalScaleFactor_!=1.0) {
456      handler_->message(CLP_BARRIER_SCALING,messages_)
457        <<"diagonal"<<diagonalScaleFactor_
458        <<CoinMessageEol;
459      diagonalNorm_*=diagonalScaleFactor_;
460    } 
461    multiplyAdd(NULL,numberTotal,0.0,diagonal_,
462                diagonalScaleFactor_);
463    int * rowsDroppedThisTime = new int [numberRows_];
464    newDropped=cholesky_->factorize(diagonal_,rowsDroppedThisTime);
465    if (newDropped) {
466      if (newDropped==-1) {
467        printf("Out of memory\n");
468        problemStatus_=4;
469        //delete all temporary regions
470        deleteWorkingData();
471        if (saveMatrix) {
472          // restore normal copy
473          delete matrix_;
474          matrix_ = saveMatrix;
475        }
476        return -1;
477      } else {
478#ifndef NDEBUG
479        //int newDropped2=cholesky_->factorize(diagonal_,rowsDroppedThisTime);
480        //assert(!newDropped2);
481#endif
482        if (newDropped<0&&0) {
483          //replace dropped
484          newDropped=-newDropped;
485          //off 1 to allow for reset all
486          newDropped--;
487          //set all bits false
488          cholesky_->resetRowsDropped();
489        } 
490      }
491    } 
492    delete [] rowsDroppedThisTime;
493    if (cholesky_->status()) {
494      std::cout << "bad cholesky?" <<std::endl;
495      abort();
496    }
497    int phase=0; // predictor, corrector , primal dual
498    double directionAccuracy=0.0;
499    bool doCorrector=true;
500    bool goodMove=true;
501    //set up for affine direction
502    setupForSolve(phase);
503    if ((modeSwitch&2)==0) {
504      directionAccuracy=findDirectionVector(phase);
505      if (directionAccuracy>worstDirectionAccuracy_) {
506        worstDirectionAccuracy_=directionAccuracy;
507      }
508      if (saveIteration>0&&directionAccuracy>1.0) {
509        handler_->message(CLP_BARRIER_EXIT2,messages_)
510          <<saveIteration
511          <<CoinMessageEol;
512        break;
513      }
514      findStepLength(phase);
515      nextGap=complementarityGap(nextNumber,nextNumberItems,1);
516      debugMove(0,actualPrimalStep_,actualDualStep_);
517      debugMove(0,1.0e-2,1.0e-2);
518    }
519    double affineGap=nextGap;
520    int bestPhase=0;
521    double bestNextGap=nextGap;
522    // ?
523    bestNextGap=CoinMax(nextGap,0.8*complementarityGap_);
524    if (quadraticObj)
525      bestNextGap=CoinMax(nextGap,0.99*complementarityGap_);
526    if (complementarityGap_>1.0e-4*numberComplementarityPairs_) {
527      //std::cout <<"predicted duality gap "<<nextGap<<std::endl;
528      double part1=nextGap/numberComplementarityPairs_;
529      part1=nextGap/numberComplementarityItems_;
530      double part2=nextGap/complementarityGap_;
531      mu_=part1*part2*part2;
532#if 0
533      double papermu =complementarityGap_/numberComplementarityPairs_;
534      double affmu = nextGap/nextNumber;
535      double sigma = pow(affmu/papermu,3);
536      printf("mu %g, papermu %g, affmu %g, sigma %g sigmamu %g\n",
537             mu_,papermu,affmu,sigma,sigma*papermu);
538#endif
539      //printf("paper mu %g\n",(nextGap*nextGap*nextGap)/(complementarityGap_*complementarityGap_*
540      //                                            (double) numberComplementarityPairs_));
541    } else {
542      double phi;
543      if (numberComplementarityPairs_<=5000) {
544        phi=pow(static_cast<double> (numberComplementarityPairs_),2.0);
545      } else {
546        phi=pow(static_cast<double> (numberComplementarityPairs_),1.5);
547        if (phi<500.0*500.0) {
548          phi=500.0*500.0;
549        } 
550      }
551      mu_=complementarityGap_/phi;
552    } 
553    //save information
554    double product=affineProduct();
555#if 0
556    //can we do corrector step?
557    double xx= complementarityGap_*(beta2-tau) +product;
558    if (xx>0.0) {
559      double saveMu = mu_;
560      double mu2=numberComplementarityPairs_;
561      mu2=xx/mu2;
562      if (mu2>mu_) {
563        //std::cout<<" could increase to "<<mu2<<std::endl;
564        //was mu2=mu2*0.25;
565        mu2=mu2*0.99;
566        if (mu2<mu_) {
567          mu_=mu2;
568          //std::cout<<" changing mu to "<<mu_<<std::endl;
569        } else {
570          //std::cout<<std::endl;
571        }
572      } else {
573        //std::cout<<" should decrease to "<<mu2<<std::endl;
574        mu_=0.5*mu2;
575        //std::cout<<" changing mu to "<<mu_<<std::endl;
576      }
577      handler_->message(CLP_BARRIER_MU,messages_)
578        <<saveMu<<mu_
579        <<CoinMessageEol;
580    } else {
581      //std::cout<<" bad by any standards"<<std::endl;
582    }
583#endif
584    if (complementarityGap_*(beta2-tau)+product-mu_*numberComplementarityPairs_<0.0&&0) {
585#ifdef SOME_DEBUG
586      printf("failed 1 product %.18g mu %.18g - %.18g < 0.0, nextGap %.18g\n",product,mu_,
587             complementarityGap_*(beta2-tau)+product-mu_*numberComplementarityPairs_,
588             nextGap);
589#endif
590      doCorrector=false;
591      if (nextGap>0.9*complementarityGap_||1) {
592        goodMove=false;
593        bestNextGap=COIN_DBL_MAX;
594      }
595      //double floatNumber = 2.0*numberComplementarityPairs_;
596      //floatNumber = 1.0*numberComplementarityItems_;
597      //mu_=nextGap/floatNumber;
598      handler_->message(CLP_BARRIER_INFO,messages_)
599        <<"no corrector step"
600        <<CoinMessageEol;
601    } else {
602      phase=1;
603    }
604    // If bad gap - try standard primal dual
605    if (nextGap>complementarityGap_*1.001)
606      goodMove=false;
607    if ((modeSwitch&2)!=0)
608      goodMove=false;
609    if (goodMove&&doCorrector) {
610      CoinMemcpyN(deltaX_,numberTotal,saveX);
611      CoinMemcpyN(deltaY_,numberRows_,saveY);
612      CoinMemcpyN(deltaZ_,numberTotal,saveZ);
613      CoinMemcpyN(deltaW_,numberTotal,saveW);
614      CoinMemcpyN(deltaSL_,numberTotal,saveSL);
615      CoinMemcpyN(deltaSU_,numberTotal,saveSU);
616#ifdef HALVE
617      double savePrimalStep = actualPrimalStep_;
618      double saveDualStep = actualDualStep_;
619      double saveMu = mu_;
620#endif
621      //set up for next step
622      setupForSolve(phase);
623      double directionAccuracy2=findDirectionVector(phase);
624      if (directionAccuracy2>worstDirectionAccuracy_) {
625        worstDirectionAccuracy_=directionAccuracy2;
626      } 
627      double testValue=1.0e2*directionAccuracy;
628      if (1.0e2*projectionTolerance_>testValue) {
629        testValue=1.0e2*projectionTolerance_;
630      } 
631      if (primalTolerance()>testValue) {
632        testValue=primalTolerance();
633      } 
634      if (maximumRHSError_>testValue) {
635        testValue=maximumRHSError_;
636      } 
637      if (directionAccuracy2>testValue&&numberIterations_>=-77) {
638        goodMove=false;
639#ifdef SOME_DEBUG
640        printf("accuracy %g phase 1 failed, test value %g\n",
641               directionAccuracy2,testValue);
642#endif
643      } 
644      if (goodMove) {
645        phase=1;
646        double norm = findStepLength(phase);
647        nextGap = complementarityGap(nextNumber,nextNumberItems,1);
648        debugMove(1,actualPrimalStep_,actualDualStep_);
649        //debugMove(1,1.0e-7,1.0e-7);
650        goodMove=checkGoodMove(true,bestNextGap,allowIncreasingGap);
651        if (norm<0)
652          goodMove=false;
653        if (!goodMove) {
654#ifdef SOME_DEBUG
655          printf("checkGoodMove failed\n");
656#endif
657        }
658      }
659#ifdef HALVE
660      int nHalve=0;
661      // relax test
662      bestNextGap=CoinMax(bestNextGap,0.9*complementarityGap_);
663      while (!goodMove) {
664        mu_=saveMu;
665        actualPrimalStep_ = savePrimalStep;
666        actualDualStep_ = saveDualStep;
667        int i;
668        //printf("halve %d\n",nHalve);
669        nHalve++;
670        const double lambda=0.5;
671        for (i=0;i<numberRows_;i++)
672          deltaY_[i] = lambda*deltaY_[i]+(1.0-lambda)*saveY[i];
673        for (i=0;i<numberTotal;i++) {
674          deltaX_[i] = lambda*deltaX_[i]+(1.0-lambda)*saveX[i];
675          deltaZ_[i] = lambda*deltaZ_[i]+(1.0-lambda)*saveZ[i];
676          deltaW_[i] = lambda*deltaW_[i]+(1.0-lambda)*saveW[i];
677          deltaSL_[i] = lambda*deltaSL_[i]+(1.0-lambda)*saveSL[i];
678          deltaSU_[i] = lambda*deltaSU_[i]+(1.0-lambda)*saveSU[i];
679        }
680   CoinMemcpyN(saveX,numberTotal,deltaX_);
681   CoinMemcpyN(saveY,numberRows_,deltaY_);
682   CoinMemcpyN(saveZ,numberTotal,deltaZ_);
683   CoinMemcpyN(saveW,numberTotal,deltaW_);
684   CoinMemcpyN(saveSL,numberTotal,deltaSL_);
685   CoinMemcpyN(saveSU,numberTotal,deltaSU_);
686        findStepLength(1);
687        nextGap = complementarityGap(nextNumber,nextNumberItems,1);
688        goodMove=checkGoodMove(true,bestNextGap,allowIncreasingGap);
689        if (nHalve>10)
690          break;
691        //assert (goodMove);
692      }
693      if (nHalve&&handler_->logLevel()>1) 
694        printf("halved %d times\n",nHalve);
695#endif
696    }
697    //bestPhase=-1;
698    //goodMove=false;
699    if (!goodMove) {
700      // Just primal dual step
701      double floatNumber;
702      floatNumber = 2.0*numberComplementarityPairs_;
703      //floatNumber = numberComplementarityItems_;
704      double saveMu=mu_; // use one from predictor corrector
705      mu_=complementarityGap_/floatNumber;
706      // If going well try small mu
707      mu_ *= sqrt((1.0-lastStep)/(1.0+10.0*lastStep));
708      double mu1=mu_;
709      double phi;
710      if (numberComplementarityPairs_<=500) {
711        phi=pow(static_cast<double> (numberComplementarityPairs_),2.0);
712      } else {
713        phi=pow(static_cast<double> (numberComplementarityPairs_),1.5);
714        if (phi<500.0*500.0) {
715          phi=500.0*500.0;
716        } 
717      }
718      mu_=complementarityGap_/phi;
719      //printf("pd mu %g, alternate %g, smallest %g\n",
720      //     mu_,mu1,smallestPrimalDualMu);
721      mu_ = sqrt(mu_*mu1);
722      mu_=mu1;
723      if ((numberIterations_&1)==0||numberIterations_<10)
724        mu_=saveMu;
725      //mu_=CoinMin(smallestPrimalDualMu*0.95,mu_);
726      smallestPrimalDualMu = mu_;
727      // Try simpler
728      floatNumber = numberComplementarityItems_;
729      mu_=0.5*complementarityGap_/floatNumber;
730      //if ((modeSwitch&2)==0) {
731      //if ((numberIterations_&1)==0)
732      //  mu_ *= 0.5;
733      //} else {
734      //mu_ *= 0.8;
735      //}
736      //set up for next step
737      setupForSolve(2);
738      findDirectionVector(2);
739      double norm=findStepLength(2);
740      // just for debug
741      bestNextGap = complementarityGap_*1.0005;
742      //bestNextGap=COIN_DBL_MAX;
743      nextGap=complementarityGap(nextNumber,nextNumberItems,2);
744      debugMove(2,actualPrimalStep_,actualDualStep_);
745      //debugMove(2,1.0e-7,1.0e-7);
746      checkGoodMove(false, bestNextGap,allowIncreasingGap);
747      if ((nextGap>0.9*complementarityGap_&&bestPhase==0&&affineGap<nextGap
748          &&(numberIterations_>80||(numberIterations_>20&&quadraticObj)))||norm<0.0) {
749        // Back to affine
750        phase=0;
751        setupForSolve(phase);
752        directionAccuracy=findDirectionVector(phase);
753        findStepLength(phase);
754        nextGap=complementarityGap(nextNumber,nextNumberItems,1);
755        bestNextGap = complementarityGap_;
756        //checkGoodMove(false, bestNextGap,allowIncreasingGap);
757      }
758    }
759    if (numberIterations_==0)
760      smallestPrimalDualMu=mu_;
761    if (!goodMove)
762      mu_=nextGap / (static_cast<double> (nextNumber)*1.1);
763    //if (quadraticObj)
764    //goodMove=true;
765    //goodMove=false; //TEMP
766    // Do centering steps
767    int numberTries=0;
768    double nextCenterGap=0.0;
769    int numberGoodTries=0;
770    double originalDualStep=actualDualStep_;
771    double originalPrimalStep=actualPrimalStep_;
772    if (actualDualStep_>0.9&&actualPrimalStep_>0.9)
773      goodMove=false; // don't bother
774    if ((modeSwitch&1)!=0)
775      goodMove=false;
776    while (goodMove&&numberTries<5) {
777      goodMove=false;
778      numberTries++;
779      CoinMemcpyN(deltaX_,numberTotal,saveX);
780      CoinMemcpyN(deltaY_,numberRows_,saveY);
781      CoinMemcpyN(deltaZ_,numberTotal,saveZ);
782      CoinMemcpyN(deltaW_,numberTotal,saveW);
783      double savePrimalStep = actualPrimalStep_;
784      double saveDualStep = actualDualStep_;
785      double saveMu = mu_;
786      setupForSolve(3);
787      findDirectionVector(3);
788      findStepLength(3);
789      debugMove(3,actualPrimalStep_,actualDualStep_);
790      //debugMove(3,1.0e-7,1.0e-7);
791      double xGap = complementarityGap(nextNumber,nextNumberItems,3);
792      // If one small then that's the one that counts
793      double checkDual=saveDualStep;
794      double checkPrimal=savePrimalStep;
795      if (checkDual>5.0*checkPrimal) {
796        checkDual=2.0*checkPrimal;
797      } else if (checkPrimal>5.0*checkDual) {
798        checkPrimal=2.0*checkDual;
799      }
800      if (actualPrimalStep_<checkPrimal||
801          actualDualStep_<checkDual||
802          (xGap>nextGap&&xGap>0.9*complementarityGap_)) {
803        //if (actualPrimalStep_<=checkPrimal||
804        //actualDualStep_<=checkDual) {
805#ifdef SOME_DEBUG
806        printf("PP rejected gap %.18g, steps %.18g %.18g, 2 gap %.18g, steps %.18g %.18g\n",xGap,
807               actualPrimalStep_,actualDualStep_,nextGap,savePrimalStep,saveDualStep);
808#endif
809        mu_=saveMu;
810        actualPrimalStep_ = savePrimalStep;
811        actualDualStep_ = saveDualStep;
812 CoinMemcpyN(saveX,numberTotal,deltaX_);
813 CoinMemcpyN(saveY,numberRows_,deltaY_);
814 CoinMemcpyN(saveZ,numberTotal,deltaZ_);
815 CoinMemcpyN(saveW,numberTotal,deltaW_);
816      } else {
817#ifdef SOME_DEBUG
818        printf("PPphase 3 gap %.18g, steps %.18g %.18g, 2 gap %.18g, steps %.18g %.18g\n",xGap,
819               actualPrimalStep_,actualDualStep_,nextGap,savePrimalStep,saveDualStep);
820#endif
821        numberGoodTries++;
822        nextCenterGap=xGap;
823        // See if big enough change
824        if (actualPrimalStep_<1.01*checkPrimal||
825            actualDualStep_<1.01*checkDual) {
826          // stop now
827        } else {
828          // carry on
829          goodMove=true;
830        }
831      }
832    }
833    if (numberGoodTries&&handler_->logLevel()>1) {
834      printf("%d centering steps moved from (gap %.18g, dual %.18g, primal %.18g) to (gap %.18g, dual %.18g, primal %.18g)\n",
835             numberGoodTries,nextGap,originalDualStep,originalPrimalStep,
836             nextCenterGap, actualDualStep_,actualPrimalStep_);
837    }
838    // save last gap
839    checkGap = complementarityGap_;
840    numberFixed=updateSolution(nextGap);
841    numberFixedTotal+=numberFixed;
842  } /* endwhile */
843  delete [] saveX;
844  delete [] saveY;
845  delete [] saveZ;
846  delete [] saveW;
847  delete [] saveSL;
848  delete [] saveSU;
849  if (savePi) {
850    std::cout<<"Restoring from iteration "<<saveIteration<<std::endl;
851    CoinMemcpyN(savePi,numberRows_,dual_);
852    CoinMemcpyN(savePrimal,numberTotal,solution_);
853    delete [] savePi;
854    delete [] savePrimal;
855  } 
856  //recompute slacks
857  // Split out solution
858  CoinZeroN(rowActivity_,numberRows_);
859  CoinMemcpyN(solution_,numberColumns_,columnActivity_);
860  matrix_->times(1.0,columnActivity_,rowActivity_);
861  //unscale objective
862  multiplyAdd(NULL,numberTotal,0.0,cost_,scaleFactor_);
863  multiplyAdd(NULL,numberRows_,0,dual_,scaleFactor_);
864  checkSolution();
865  CoinMemcpyN(reducedCost_,numberColumns_,dj_);
866  // If quadratic use last solution
867  // Restore quadratic objective if necessary
868  if (saveObjective) {
869    delete objective_;
870    objective_=saveObjective;
871    objectiveValue_=0.5*(primalObjective_+dualObjective_);
872  }
873  handler_->message(CLP_BARRIER_END,messages_)
874    <<sumPrimalInfeasibilities_
875    <<sumDualInfeasibilities_
876    <<complementarityGap_
877    <<objectiveValue()
878    <<CoinMessageEol;
879  //#ifdef SOME_DEBUG
880  if (handler_->logLevel()>1) 
881    printf("ENDRUN status %d after %d iterations\n",problemStatus_,numberIterations_);
882  //#endif
883  //std::cout<<"Absolute primal infeasibility at end "<<sumPrimalInfeasibilities_<<std::endl;
884  //std::cout<<"Absolute dual infeasibility at end "<<sumDualInfeasibilities_<<std::endl;
885  //std::cout<<"Absolute complementarity at end "<<complementarityGap_<<std::endl;
886  //std::cout<<"Primal objective "<<objectiveValue()<<std::endl;
887  //std::cout<<"maximum complementarity "<<worstComplementarity_<<std::endl;
888  //delete all temporary regions
889  deleteWorkingData();
890  if (saveMatrix) {
891    // restore normal copy
892    delete matrix_;
893    matrix_ = saveMatrix;
894  }
895  return problemStatus_;
896}
897// findStepLength.
898//phase  - 0 predictor
899//         1 corrector
900//         2 primal dual
901double ClpPredictorCorrector::findStepLength( int phase)
902{
903  double directionNorm=0.0;
904  double maximumPrimalStep=COIN_DBL_MAX;
905  double maximumDualStep=COIN_DBL_MAX;
906  int numberTotal = numberRows_+numberColumns_;
907  double tolerance = 1.0e-12;
908  int chosenPrimalSequence=-1;
909  int chosenDualSequence=-1;
910  bool lowPrimal=false;
911  bool lowDual=false;
912  // If done many iterations then allow to hit boundary
913  double hitTolerance;
914  //printf("objective norm %g\n",objectiveNorm_);
915  if (numberIterations_<80||!gonePrimalFeasible_)
916    hitTolerance = COIN_DBL_MAX;
917  else
918    hitTolerance = CoinMax(1.0e3,1.0e-3*objectiveNorm_);
919  int iColumn;
920  //printf("dual value %g\n",dual_[0]);
921  //printf("     X     dX      lS     dlS     uS     dUs    dj    Z dZ     t   dT\n");
922  for (iColumn=0;iColumn<numberTotal;iColumn++) {
923    if (!flagged(iColumn)) {
924      double directionElement=deltaX_[iColumn];
925      if (directionNorm<fabs(directionElement)) {
926        directionNorm=fabs(directionElement);
927      }
928      if (0) 
929      printf("%d %g %g %g %g %g %g %g %g %g %g %g\n",
930             iColumn,solution_[iColumn],deltaX_[iColumn],
931             lowerSlack_[iColumn],deltaSL_[iColumn],
932             upperSlack_[iColumn],deltaSU_[iColumn],
933             dj_[iColumn],
934             zVec_[iColumn],deltaZ_[iColumn],
935             wVec_[iColumn],deltaW_[iColumn]);
936      if (lowerBound(iColumn)) {
937        double delta = - deltaSL_[iColumn];
938        double z1 = deltaZ_[iColumn];
939        double newZ = zVec_[iColumn]+z1;
940        if (zVec_[iColumn]>tolerance) {
941          if (zVec_[iColumn]<-z1*maximumDualStep) {
942            maximumDualStep=-zVec_[iColumn]/z1;
943            chosenDualSequence=iColumn;
944            lowDual=true;
945          } 
946        } 
947        if (lowerSlack_[iColumn]<maximumPrimalStep*delta) {
948          double newStep=lowerSlack_[iColumn]/delta;
949          if (newStep>0.2||newZ<hitTolerance||delta>1.0e3||delta<=1.0e-6||dj_[iColumn]<hitTolerance) {
950            maximumPrimalStep = newStep;
951            chosenPrimalSequence=iColumn;
952            lowPrimal=true;
953          } else {
954            //printf("small %d delta %g newZ %g step %g\n",iColumn,delta,newZ,newStep);
955          }
956        } 
957      }
958      if (upperBound(iColumn)) {
959        double delta = - deltaSU_[iColumn];;
960        double w1 = deltaW_[iColumn];
961        double newT = wVec_[iColumn]+w1;
962        if (wVec_[iColumn]>tolerance) {
963          if (wVec_[iColumn]<-w1*maximumDualStep) {
964            maximumDualStep=-wVec_[iColumn]/w1;
965            chosenDualSequence=iColumn;
966            lowDual=false;
967          } 
968        } 
969        if (upperSlack_[iColumn]<maximumPrimalStep*delta) {
970          double newStep=upperSlack_[iColumn]/delta;
971          if (newStep>0.2||newT<hitTolerance||delta>1.0e3||delta<=1.0e-6||dj_[iColumn]>-hitTolerance) {
972            maximumPrimalStep = newStep;
973            chosenPrimalSequence=iColumn;
974            lowPrimal=false;
975          } else {
976            //printf("small %d delta %g newT %g step %g\n",iColumn,delta,newT,newStep);
977          }
978        } 
979      } 
980    } 
981  }
982#ifdef SOME_DEBUG
983  printf("new step - phase %d, norm %.18g, dual step %.18g, primal step %.18g\n",
984         phase,directionNorm,maximumDualStep,maximumPrimalStep);
985  if (lowDual) 
986    printf("ld %d %g %g => %g (dj %g,sol %g) ",
987           chosenDualSequence,zVec_[chosenDualSequence],
988           deltaZ_[chosenDualSequence],zVec_[chosenDualSequence]+
989           maximumDualStep*deltaZ_[chosenDualSequence],dj_[chosenDualSequence],
990           solution_[chosenDualSequence]);
991  else
992    printf("ud %d %g %g => %g (dj %g,sol %g) ",
993           chosenDualSequence,wVec_[chosenDualSequence],
994           deltaW_[chosenDualSequence],wVec_[chosenDualSequence]+
995           maximumDualStep*deltaW_[chosenDualSequence],dj_[chosenDualSequence],
996           solution_[chosenDualSequence]);
997  if (lowPrimal) 
998    printf("lp %d %g %g => %g (dj %g,sol %g)\n",
999           chosenPrimalSequence,lowerSlack_[chosenPrimalSequence],
1000           deltaSL_[chosenPrimalSequence],lowerSlack_[chosenPrimalSequence]+
1001           maximumPrimalStep*deltaSL_[chosenPrimalSequence],
1002           dj_[chosenPrimalSequence],solution_[chosenPrimalSequence]);
1003  else
1004    printf("up %d %g %g => %g (dj %g,sol %g)\n",
1005           chosenPrimalSequence,upperSlack_[chosenPrimalSequence],
1006           deltaSU_[chosenPrimalSequence],upperSlack_[chosenPrimalSequence]+
1007           maximumPrimalStep*deltaSU_[chosenPrimalSequence],
1008           dj_[chosenPrimalSequence],solution_[chosenPrimalSequence]);
1009#endif
1010  actualPrimalStep_=stepLength_*maximumPrimalStep;
1011  if (phase>=0&&actualPrimalStep_>1.0) {
1012    actualPrimalStep_=1.0;
1013  } 
1014  actualDualStep_=stepLength_*maximumDualStep;
1015  if (phase>=0&&actualDualStep_>1.0) {
1016    actualDualStep_=1.0;
1017  } 
1018  // See if quadratic objective
1019#ifndef NO_RTTI
1020  ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
1021#else
1022  ClpQuadraticObjective * quadraticObj = NULL;
1023  if (objective_->type()==2)
1024    quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
1025#endif
1026  if (quadraticObj) {
1027    // Use smaller unless very small
1028    double smallerStep=CoinMin(actualDualStep_,actualPrimalStep_);
1029    if (smallerStep>0.0001) {
1030      actualDualStep_=smallerStep;
1031      actualPrimalStep_=smallerStep;
1032    }
1033  }
1034#define OFFQ
1035#ifndef OFFQ
1036  if (quadraticObj) {
1037    // Don't bother if phase 0 or 3 or large gap
1038    //if ((phase==1||phase==2||phase==0)&&maximumDualError_>0.1*complementarityGap_
1039    //&&smallerStep>0.001) {
1040    if ((phase==1||phase==2||phase==0||phase==3)) {
1041      // minimize complementarity + norm*dual inf ? primal inf
1042      // at first - just check better - if not
1043      // Complementarity gap will be a*change*change + b*change +c
1044      double a=0.0;
1045      double b=0.0;
1046      double c=0.0;
1047      /* SQUARE of dual infeasibility will be:
1048         square of dj - ......
1049      */
1050      double aq=0.0;
1051      double bq=0.0;
1052      double cq=0.0;
1053      double gamma2 = gamma_*gamma_; // gamma*gamma will be added to diagonal
1054      double * linearDjChange = new double[numberTotal];
1055      CoinZeroN(linearDjChange,numberColumns_);
1056      multiplyAdd(deltaY_,numberRows_,1.0,linearDjChange+numberColumns_,0.0);
1057      matrix_->transposeTimes(-1.0,deltaY_,linearDjChange);
1058      CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
1059      const int * columnQuadratic = quadratic->getIndices();
1060      const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
1061      const int * columnQuadraticLength = quadratic->getVectorLengths();
1062      double * quadraticElement = quadratic->getMutableElements();
1063      for (iColumn=0;iColumn<numberTotal;iColumn++) {
1064        double oldPrimal = solution_[iColumn];
1065        if (!flagged(iColumn)) {
1066          if (lowerBound(iColumn)) {
1067            double change =oldPrimal+deltaX_[iColumn]-lowerSlack_[iColumn]-lower_[iColumn];
1068            c += lowerSlack_[iColumn]*zVec_[iColumn];
1069            b += lowerSlack_[iColumn]*deltaZ_[iColumn]+zVec_[iColumn]*change;
1070            a += deltaZ_[iColumn]*change;
1071          }
1072          if (upperBound(iColumn)) {
1073            double change =upper_[iColumn]-oldPrimal-deltaX_[iColumn]-upperSlack_[iColumn];
1074            c += upperSlack_[iColumn]*wVec_[iColumn];
1075            b += upperSlack_[iColumn]*deltaW_[iColumn]+wVec_[iColumn]*change;
1076            a += deltaW_[iColumn]*change;
1077          } 
1078          // new djs are dj_ + change*value
1079          double djChange = linearDjChange[iColumn];
1080          if (iColumn<numberColumns_) {
1081            for (CoinBigIndex j=columnQuadraticStart[iColumn];
1082                 j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
1083              int jColumn = columnQuadratic[j];
1084              double changeJ = deltaX_[jColumn];
1085              double elementValue = quadraticElement[j];
1086              djChange += changeJ*elementValue;
1087            }
1088          }
1089          double gammaTerm = gamma2;
1090          if (primalR_) {
1091            gammaTerm += primalR_[iColumn];
1092          }
1093          djChange += gammaTerm;
1094          // and dual infeasibility
1095          double oldInf = dj_[iColumn]-zVec_[iColumn]+wVec_[iColumn]+
1096            gammaTerm*solution_[iColumn];
1097          double changeInf = djChange-deltaZ_[iColumn]+deltaW_[iColumn];
1098          cq += oldInf*oldInf;
1099          bq += 2.0*oldInf*changeInf;
1100          aq += changeInf*changeInf;
1101        } else {
1102          // fixed
1103          if (lowerBound(iColumn)) {
1104            c += lowerSlack_[iColumn]*zVec_[iColumn];
1105          }
1106          if (upperBound(iColumn)) {
1107            c += upperSlack_[iColumn]*wVec_[iColumn];
1108          } 
1109          // new djs are dj_ + change*value
1110          double djChange = linearDjChange[iColumn];
1111          if (iColumn<numberColumns_) {
1112            for (CoinBigIndex j=columnQuadraticStart[iColumn];
1113                 j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
1114              int jColumn = columnQuadratic[j];
1115              double changeJ = deltaX_[jColumn];
1116              double elementValue = quadraticElement[j];
1117              djChange += changeJ*elementValue;
1118            }
1119          }
1120          double gammaTerm = gamma2;
1121          if (primalR_) {
1122            gammaTerm += primalR_[iColumn];
1123          }
1124          djChange += gammaTerm;
1125          // and dual infeasibility
1126          double oldInf = dj_[iColumn]-zVec_[iColumn]+wVec_[iColumn]+
1127            gammaTerm*solution_[iColumn];
1128          double changeInf = djChange-deltaZ_[iColumn]+deltaW_[iColumn];
1129          cq += oldInf*oldInf;
1130          bq += 2.0*oldInf*changeInf;
1131          aq += changeInf*changeInf;
1132        }
1133      }
1134      delete [] linearDjChange;
1135      // ? We want to minimize complementarityGap + solutionNorm_*square of inf ??
1136      // maybe use inf and do line search
1137      // To check see if matches at current step
1138      double step=actualPrimalStep_;
1139      //Current gap + solutionNorm_ * sqrt (sum square inf)
1140      double multiplier = solutionNorm_;
1141      multiplier *= 0.01;
1142      multiplier=1.0;
1143      double currentInf =  multiplier*sqrt(cq);
1144      double nextInf =  multiplier*sqrt(CoinMax(cq+step*bq+step*step*aq,0.0));
1145      double allowedIncrease=1.4;
1146#ifdef SOME_DEBUG
1147      printf("lin %g %g %g -> %g\n",a,b,c,
1148             c+b*step+a*step*step);
1149      printf("quad %g %g %g -> %g\n",aq,bq,cq,
1150             cq+bq*step+aq*step*step);
1151      debugMove(7,step,step);
1152      printf ("current dualInf %g, with step of %g is %g\n",
1153              currentInf,step,nextInf);
1154#endif
1155      if (b>-1.0e-6) {
1156        if (phase!=0)
1157          directionNorm=-1.0;
1158      }
1159      if ((phase==1||phase==2||phase==0||phase==3)&&nextInf>0.1*complementarityGap_&&
1160          nextInf>currentInf*allowedIncrease) {
1161        //cq = CoinMax(cq,10.0);
1162        // convert to (x+q)*(x+q) = w
1163        double q = bq/(1.0*aq);
1164        double w = CoinMax(q*q + (cq/aq)*(allowedIncrease-1.0),0.0);
1165        w = sqrt(w);
1166        double stepX = w-q;
1167        step=stepX;
1168        nextInf = 
1169          multiplier*sqrt(CoinMax(cq+step*bq+step*step*aq,0.0));
1170#ifdef SOME_DEBUG
1171        printf ("with step of %g dualInf is %g\n",
1172                step,nextInf);
1173#endif
1174        actualDualStep_=CoinMin(step,actualDualStep_);
1175        actualPrimalStep_=CoinMin(step,actualPrimalStep_);
1176      }
1177    }
1178  } else {
1179    // probably pointless as linear
1180    // minimize complementarity
1181    // Complementarity gap will be a*change*change + b*change +c
1182    double a=0.0;
1183    double b=0.0;
1184    double c=0.0;
1185    for (iColumn=0;iColumn<numberTotal;iColumn++) {
1186      double oldPrimal = solution_[iColumn];
1187      if (!flagged(iColumn)) {
1188        if (lowerBound(iColumn)) {
1189          double change =oldPrimal+deltaX_[iColumn]-lowerSlack_[iColumn]-lower_[iColumn];
1190          c += lowerSlack_[iColumn]*zVec_[iColumn];
1191          b += lowerSlack_[iColumn]*deltaZ_[iColumn]+zVec_[iColumn]*change;
1192          a += deltaZ_[iColumn]*change;
1193        }
1194        if (upperBound(iColumn)) {
1195          double change =upper_[iColumn]-oldPrimal-deltaX_[iColumn]-upperSlack_[iColumn];
1196          c += upperSlack_[iColumn]*wVec_[iColumn];
1197          b += upperSlack_[iColumn]*deltaW_[iColumn]+wVec_[iColumn]*change;
1198          a += deltaW_[iColumn]*change;
1199        } 
1200      } else {
1201        // fixed
1202        if (lowerBound(iColumn)) {
1203          c += lowerSlack_[iColumn]*zVec_[iColumn];
1204        }
1205        if (upperBound(iColumn)) {
1206          c += upperSlack_[iColumn]*wVec_[iColumn];
1207        } 
1208      }
1209    }
1210    // ? We want to minimize complementarityGap;
1211    // maybe use inf and do line search
1212    // To check see if matches at current step
1213    double step=CoinMin(actualPrimalStep_,actualDualStep_);
1214    double next = c+b*step+a*step*step;
1215#ifdef SOME_DEBUG
1216    printf("lin %g %g %g -> %g\n",a,b,c,
1217           c+b*step+a*step*step);
1218    debugMove(7,step,step);
1219#endif
1220    if (b>-1.0e-6) {
1221      if (phase==0) {
1222#ifdef SOME_DEBUG
1223        printf("*** odd phase 0 direction\n");
1224#endif
1225      } else {
1226        directionNorm=-1.0;
1227      }
1228    }
1229    // and with ratio
1230    a=0.0;
1231    b=0.0;
1232    double ratio = actualDualStep_/actualPrimalStep_;
1233    for (iColumn=0;iColumn<numberTotal;iColumn++) {
1234      double oldPrimal = solution_[iColumn];
1235      if (!flagged(iColumn)) {
1236        if (lowerBound(iColumn)) {
1237          double change =oldPrimal+deltaX_[iColumn]-lowerSlack_[iColumn]-lower_[iColumn];
1238          b += lowerSlack_[iColumn]*deltaZ_[iColumn]*ratio+zVec_[iColumn]*change;
1239          a += deltaZ_[iColumn]*change*ratio;
1240        }
1241        if (upperBound(iColumn)) {
1242          double change =upper_[iColumn]-oldPrimal-deltaX_[iColumn]-upperSlack_[iColumn];
1243          b += upperSlack_[iColumn]*deltaW_[iColumn]*ratio+wVec_[iColumn]*change;
1244          a += deltaW_[iColumn]*change*ratio;
1245        } 
1246      }
1247    }
1248    // ? We want to minimize complementarityGap;
1249    // maybe use inf and do line search
1250    // To check see if matches at current step
1251    step=actualPrimalStep_;
1252    double next2 = c+b*step+a*step*step;
1253    if (next2>next) {
1254      actualPrimalStep_=CoinMin(actualPrimalStep_,actualDualStep_);
1255      actualDualStep_=actualPrimalStep_;
1256    }
1257#ifdef SOME_DEBUG
1258    printf("linb %g %g %g -> %g\n",a,b,c,
1259           c+b*step+a*step*step);
1260    debugMove(7,actualPrimalStep_,actualDualStep_);
1261#endif
1262    if (b>-1.0e-6) {
1263      if (phase==0) {
1264#ifdef SOME_DEBUG
1265        printf("*** odd phase 0 direction\n");
1266#endif
1267      } else {
1268        directionNorm=-1.0;
1269      }
1270    }
1271  }
1272#else
1273  //actualPrimalStep_ =0.5*actualDualStep_;
1274#endif
1275#ifdef FULL_DEBUG
1276  if (phase==3){
1277    double minBeta = 0.1*mu_;
1278    double maxBeta = 10.0*mu_;
1279    for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
1280      if (!flagged(iColumn)) {
1281        if (lowerBound(iColumn)) {
1282          double change = -rhsL_[iColumn] + deltaX_[iColumn];
1283          double dualValue=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn];
1284          double primalValue=lowerSlack_[iColumn]+actualPrimalStep_*change;
1285          double gapProduct=dualValue*primalValue;
1286          if (delta2Z_[iColumn]<minBeta||delta2Z_[iColumn]>maxBeta)
1287            printf("3lower %d primal %g, dual %g, gap %g, old gap %g\n",
1288                   iColumn,primalValue,dualValue,gapProduct,delta2Z_[iColumn]);
1289        } 
1290        if (upperBound(iColumn)) {
1291          double change = rhsU_[iColumn]-deltaX_[iColumn];
1292          double dualValue=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn];
1293          double primalValue=upperSlack_[iColumn]+actualPrimalStep_*change;
1294          double gapProduct=dualValue*primalValue;
1295          if (delta2W_[iColumn]<minBeta||delta2W_[iColumn]>maxBeta)
1296            printf("3upper %d primal %g, dual %g, gap %g, old gap %g\n",
1297                 iColumn,primalValue,dualValue,gapProduct,delta2W_[iColumn]);
1298        } 
1299      } 
1300    }
1301  }
1302#endif
1303#ifdef SOME_DEBUG_not
1304  {
1305    double largestL=0.0;
1306    double smallestL=COIN_DBL_MAX;
1307    double largestU=0.0;
1308    double smallestU=COIN_DBL_MAX;
1309    double sumL=0.0;
1310    double sumU=0.0;
1311    int nL=0;
1312    int nU=0;
1313    for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
1314      if (!flagged(iColumn)) {
1315        if (lowerBound(iColumn)) {
1316          double change = -rhsL_[iColumn] + deltaX_[iColumn];
1317          double dualValue=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn];
1318          double primalValue=lowerSlack_[iColumn]+actualPrimalStep_*change;
1319          double gapProduct=dualValue*primalValue;
1320          largestL = CoinMax(largestL,gapProduct);
1321          smallestL = CoinMin(smallestL,gapProduct);
1322          nL++;
1323          sumL += gapProduct;
1324        } 
1325        if (upperBound(iColumn)) {
1326          double change = rhsU_[iColumn]-deltaX_[iColumn];
1327          double dualValue=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn];
1328          double primalValue=upperSlack_[iColumn]+actualPrimalStep_*change;
1329          double gapProduct=dualValue*primalValue;
1330          largestU = CoinMax(largestU,gapProduct);
1331          smallestU = CoinMin(smallestU,gapProduct);
1332          nU++;
1333          sumU += gapProduct;
1334        } 
1335      } 
1336    }
1337    double mu = (sumL+sumU)/(static_cast<double> (nL+nU));
1338
1339    double minBeta = 0.1*mu;
1340    double maxBeta = 10.0*mu;
1341    int nBL=0;
1342    int nAL=0;
1343    int nBU=0;
1344    int nAU=0;
1345    for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
1346      if (!flagged(iColumn)) {
1347        if (lowerBound(iColumn)) {
1348          double change = -rhsL_[iColumn] + deltaX_[iColumn];
1349          double dualValue=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn];
1350          double primalValue=lowerSlack_[iColumn]+actualPrimalStep_*change;
1351          double gapProduct=dualValue*primalValue;
1352          if (gapProduct<minBeta)
1353            nBL++;
1354          else if (gapProduct>maxBeta)
1355            nAL++;
1356          //if (gapProduct<0.1*minBeta)
1357          //printf("Lsmall one %d dual %g primal %g\n",iColumn,
1358          //   dualValue,primalValue);
1359        } 
1360        if (upperBound(iColumn)) {
1361          double change = rhsU_[iColumn]-deltaX_[iColumn];
1362          double dualValue=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn];
1363          double primalValue=upperSlack_[iColumn]+actualPrimalStep_*change;
1364          double gapProduct=dualValue*primalValue;
1365          if (gapProduct<minBeta)
1366            nBU++;
1367          else if (gapProduct>maxBeta)
1368            nAU++;
1369          //if (gapProduct<0.1*minBeta)
1370          //printf("Usmall one %d dual %g primal %g\n",iColumn,
1371          //   dualValue,primalValue);
1372        } 
1373      } 
1374    }
1375    printf("phase %d new mu %.18g new gap %.18g\n",phase,mu,sumL+sumU);
1376    printf("          %d lower, smallest %.18g, %d below - largest %.18g, %d above\n",
1377           nL,smallestL,nBL,largestL,nAL);
1378    printf("          %d upper, smallest %.18g, %d below - largest %.18g, %d above\n",
1379           nU,smallestU,nBU,largestU,nAU);
1380  }
1381#endif
1382  return directionNorm;
1383}
1384/* Does solve. region1 is for deltaX (columns+rows), region2 for deltaPi (rows) */
1385void 
1386ClpPredictorCorrector::solveSystem(double * region1, double * region2,
1387                                   const double * region1In, const double * region2In,
1388                                   const double * saveRegion1, const double * saveRegion2,
1389                                   bool gentleRefine)
1390{
1391  int iRow;
1392  int numberTotal = numberRows_+numberColumns_;
1393  if (region2In) {
1394    // normal
1395    for (iRow=0;iRow<numberRows_;iRow++)
1396      region2[iRow] = region2In[iRow];
1397  } else {
1398    // initial solution - (diagonal is 1 or 0)
1399    CoinZeroN(region2,numberRows_);
1400  }
1401  int iColumn;
1402  if (cholesky_->type()<20) {
1403    // not KKT
1404    for (iColumn=0;iColumn<numberTotal;iColumn++)
1405      region1[iColumn] = region1In[iColumn]*diagonal_[iColumn];
1406    multiplyAdd(region1+numberColumns_,numberRows_,-1.0,region2,1.0);
1407    matrix_->times(1.0,region1,region2);
1408    double maximumRHS = maximumAbsElement(region2,numberRows_);
1409    double scale=1.0;
1410    double unscale=1.0;
1411    if (maximumRHS>1.0e-30) {
1412      if (maximumRHS<=0.5) {
1413        double factor=2.0;
1414        while (maximumRHS<=0.5) {
1415          maximumRHS*=factor;
1416          scale*=factor;
1417        } /* endwhile */
1418      } else if (maximumRHS>=2.0&&maximumRHS<=COIN_DBL_MAX) {
1419        double factor=0.5;
1420        while (maximumRHS>=2.0) {
1421          maximumRHS*=factor;
1422          scale*=factor;
1423        } /* endwhile */
1424      } 
1425      unscale=diagonalScaleFactor_/scale;
1426    } else {
1427      //effectively zero
1428      scale=0.0;
1429      unscale=0.0;
1430    } 
1431    multiplyAdd(NULL,numberRows_,0.0,region2,scale);
1432    cholesky_->solve(region2);
1433    multiplyAdd(NULL,numberRows_,0.0,region2,unscale);
1434    multiplyAdd(region2,numberRows_,-1.0,region1+numberColumns_,0.0);
1435    CoinZeroN(region1,numberColumns_);
1436    matrix_->transposeTimes(1.0,region2,region1);
1437    for (iColumn=0;iColumn<numberTotal;iColumn++)
1438      region1[iColumn] = (region1[iColumn]-region1In[iColumn])*diagonal_[iColumn];
1439  } else {
1440    for (iColumn=0;iColumn<numberTotal;iColumn++)
1441      region1[iColumn] = region1In[iColumn];
1442    cholesky_->solveKKT(region1,region2,diagonal_,diagonalScaleFactor_);
1443  }
1444  if (saveRegion2) {
1445    //refine?
1446    double scaleX=1.0;
1447    if (gentleRefine) 
1448      scaleX=0.8;
1449    multiplyAdd(saveRegion2,numberRows_,1.0,region2,scaleX);
1450    assert (saveRegion1);
1451    multiplyAdd(saveRegion1,numberTotal,1.0,region1,scaleX);
1452  } 
1453}
1454// findDirectionVector.
1455double ClpPredictorCorrector::findDirectionVector(const int phase)
1456{
1457  double projectionTolerance=projectionTolerance_;
1458  //temporary
1459  //projectionTolerance=1.0e-15;
1460  double errorCheck=0.9*maximumRHSError_/solutionNorm_;
1461  if (errorCheck>primalTolerance()) {
1462    if (errorCheck<projectionTolerance) {
1463      projectionTolerance=errorCheck;
1464    } 
1465  } else {
1466    if (primalTolerance()<projectionTolerance) {
1467      projectionTolerance=primalTolerance();
1468    } 
1469  } 
1470  double * newError = new double [numberRows_];
1471  int numberTotal = numberRows_+numberColumns_;
1472  //if flagged then entries zero so can do
1473  // For KKT separate out
1474  double * region1Save=NULL;//for refinement
1475  int iColumn;
1476  if (cholesky_->type()<20) {
1477    int iColumn;
1478    for (iColumn=0;iColumn<numberTotal;iColumn++)
1479      deltaX_[iColumn] = workArray_[iColumn] - solution_[iColumn];
1480    multiplyAdd(deltaX_+numberColumns_,numberRows_,-1.0,deltaY_,0.0);
1481    matrix_->times(1.0,deltaX_,deltaY_);
1482  } else {
1483    // regions in will be workArray and newError
1484    // regions out deltaX_ and deltaY_
1485    multiplyAdd(solution_+numberColumns_,numberRows_,1.0,newError,0.0);
1486    matrix_->times(-1.0,solution_,newError);
1487    // This is inefficient but just for now get values which will be in deltay
1488    int iColumn;
1489    for (iColumn=0;iColumn<numberTotal;iColumn++)
1490      deltaX_[iColumn] = workArray_[iColumn] - solution_[iColumn];
1491    multiplyAdd(deltaX_+numberColumns_,numberRows_,-1.0,deltaY_,0.0);
1492    matrix_->times(1.0,deltaX_,deltaY_);
1493  }
1494  bool goodSolve=false;
1495  double * regionSave=NULL;//for refinement
1496  int numberTries=0;
1497  double relativeError=COIN_DBL_MAX;
1498  double tryError=1.0e31;
1499  while (!goodSolve&&numberTries<30) {
1500    double lastError=relativeError;
1501    goodSolve=true;
1502    double maximumRHS;
1503    double saveMaximum;
1504    maximumRHS = CoinMax(maximumAbsElement(deltaY_,numberRows_),1.0e-12);
1505    saveMaximum = maximumRHS;
1506    if (cholesky_->type()<20) {
1507      // no kkt
1508      double scale=1.0;
1509      double unscale=1.0;
1510      if (maximumRHS>1.0e-30) {
1511        if (maximumRHS<=0.5) {
1512          double factor=2.0;
1513          while (maximumRHS<=0.5) {
1514            maximumRHS*=factor;
1515            scale*=factor;
1516          } /* endwhile */
1517        } else if (maximumRHS>=2.0&&maximumRHS<=COIN_DBL_MAX) {
1518          double factor=0.5;
1519          while (maximumRHS>=2.0) {
1520            maximumRHS*=factor;
1521            scale*=factor;
1522          } /* endwhile */
1523        } 
1524        unscale=diagonalScaleFactor_/scale;
1525      } else {
1526        //effectively zero
1527        scale=0.0;
1528        unscale=0.0;
1529      } 
1530      //printf("--putting scales to 1.0\n");
1531      //scale=1.0;
1532      //unscale=1.0;
1533      multiplyAdd(NULL,numberRows_,0.0,deltaY_,scale);
1534      cholesky_->solve(deltaY_);
1535      multiplyAdd(NULL,numberRows_,0.0,deltaY_,unscale);
1536#if 0
1537      {
1538        printf("deltay\n");
1539        for (int i=0;i<numberRows_;i++)
1540          printf("%d %.18g\n",i,deltaY_[i]);
1541      }
1542      exit(66);
1543#endif
1544      if (numberTries) {
1545        //refine?
1546        double scaleX=1.0;
1547        if (lastError>1.0e-5) 
1548          scaleX=0.8;
1549        multiplyAdd(regionSave,numberRows_,1.0,deltaY_,scaleX);
1550      } 
1551      //CoinZeroN(newError,numberRows_);
1552      multiplyAdd(deltaY_,numberRows_,-1.0,deltaX_+numberColumns_,0.0);
1553      CoinZeroN(deltaX_,numberColumns_);
1554      matrix_->transposeTimes(1.0,deltaY_,deltaX_);
1555      //if flagged then entries zero so can do
1556      for (iColumn=0;iColumn<numberTotal;iColumn++)
1557        deltaX_[iColumn] = deltaX_[iColumn]*diagonal_[iColumn]
1558          -workArray_[iColumn];
1559    } else {
1560      // KKT
1561      solveSystem(deltaX_, deltaY_,
1562                  workArray_,newError,region1Save,regionSave,lastError>1.0e-5);
1563    }
1564    multiplyAdd(deltaX_+numberColumns_,numberRows_,-1.0,newError,0.0);
1565    matrix_->times(1.0,deltaX_,newError);
1566    numberTries++;
1567   
1568    //now add in old Ax - doing extra checking
1569    double maximumRHSError=0.0;
1570    double maximumRHSChange=0.0;
1571    int iRow;
1572    char * dropped = cholesky_->rowsDropped();
1573    for (iRow=0;iRow<numberRows_;iRow++) {
1574      if (!dropped[iRow]) {
1575        double newValue=newError[iRow];
1576        double oldValue=errorRegion_[iRow];
1577        //severity of errors depend on signs
1578        //**later                                                             */
1579        if (fabs(newValue)>maximumRHSChange) {
1580          maximumRHSChange=fabs(newValue);
1581        } 
1582        double result=newValue+oldValue;
1583        if (fabs(result)>maximumRHSError) {
1584          maximumRHSError=fabs(result);
1585        } 
1586        newError[iRow]=result;
1587      } else {
1588        double newValue=newError[iRow];
1589        double oldValue=errorRegion_[iRow];
1590        if (fabs(newValue)>maximumRHSChange) {
1591          maximumRHSChange=fabs(newValue);
1592        } 
1593        double result=newValue+oldValue;
1594        newError[iRow]=result;
1595        //newError[iRow]=0.0;
1596        //assert(deltaY_[iRow]==0.0);
1597        deltaY_[iRow]=0.0;
1598      } 
1599    } 
1600    relativeError = maximumRHSError/solutionNorm_;
1601    relativeError = maximumRHSError/saveMaximum;
1602    if (relativeError>tryError) 
1603      relativeError=tryError;
1604    if (relativeError<lastError) {
1605      maximumRHSChange_= maximumRHSChange;
1606      if (relativeError>1.0e-9
1607          ||numberTries>1) {
1608        handler_->message(CLP_BARRIER_ACCURACY,messages_)
1609          <<phase<<numberTries<<relativeError
1610          <<CoinMessageEol;
1611      } 
1612      if (relativeError>projectionTolerance&&numberTries<=3) {
1613        //try and refine
1614        goodSolve=false;
1615      } 
1616      //*** extra test here
1617      if (!goodSolve) {
1618        if (!regionSave) {
1619          regionSave = new double [numberRows_];
1620          if (cholesky_->type()>=20) 
1621            region1Save = new double [numberTotal];
1622        } 
1623        CoinMemcpyN(deltaY_,numberRows_,regionSave);
1624        if (cholesky_->type()<20) {
1625          // not KKT
1626          multiplyAdd(newError,numberRows_,-1.0,deltaY_,0.0);
1627        } else {
1628          // KKT
1629          CoinMemcpyN(deltaX_,numberTotal,region1Save);
1630          // and back to input region
1631          CoinMemcpyN(deltaY_,numberRows_,newError);
1632        }
1633      } 
1634    } else {
1635      //std::cout <<" worse residual = "<<relativeError;
1636      //bring back previous
1637      relativeError=lastError;
1638      if (regionSave) {
1639        CoinMemcpyN(regionSave,numberRows_,deltaY_);
1640        if (cholesky_->type()<20) {
1641          // not KKT
1642          multiplyAdd(deltaY_,numberRows_,-1.0,deltaX_+numberColumns_,0.0);
1643          CoinZeroN(deltaX_,numberColumns_);
1644          matrix_->transposeTimes(1.0,deltaY_,deltaX_);
1645          //if flagged then entries zero so can do
1646          for (iColumn=0;iColumn<numberTotal;iColumn++)
1647            deltaX_[iColumn] = deltaX_[iColumn]*diagonal_[iColumn]
1648              -workArray_[iColumn];
1649        } else {
1650          // KKT
1651          CoinMemcpyN(region1Save,numberTotal,deltaX_);
1652        }
1653      } else {
1654        // disaster
1655        CoinFillN(deltaX_,numberTotal,1.0);
1656        CoinFillN(deltaY_,numberRows_,1.0);
1657        printf("bad cholesky\n");
1658      }
1659    }
1660  } /* endwhile */
1661  delete [] regionSave;
1662  delete [] region1Save;
1663  delete [] newError;
1664  // now rest
1665  double extra=eExtra;
1666  //multiplyAdd(deltaY_,numberRows_,1.0,deltaW_+numberColumns_,0.0);
1667  //CoinZeroN(deltaW_,numberColumns_);
1668  //matrix_->transposeTimes(-1.0,deltaY_,deltaW_);
1669 
1670  for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
1671    deltaSU_[iColumn]=0.0;
1672    deltaSL_[iColumn]=0.0;
1673    deltaZ_[iColumn]=0.0;
1674    double dd=deltaW_[iColumn];
1675    deltaW_[iColumn]=0.0;
1676    if (!flagged(iColumn)) {
1677      double deltaX = deltaX_[iColumn];
1678      if (lowerBound(iColumn)) {
1679        double zValue = rhsZ_[iColumn];
1680        double gHat = zValue + zVec_[iColumn]*rhsL_[iColumn];
1681        double slack = lowerSlack_[iColumn]+extra;
1682        deltaSL_[iColumn] = -rhsL_[iColumn]+deltaX;
1683        deltaZ_[iColumn]=(gHat-zVec_[iColumn]*deltaX)/slack;
1684      } 
1685      if (upperBound(iColumn)) {
1686        double wValue = rhsW_[iColumn];
1687        double hHat = wValue - wVec_[iColumn]*rhsU_[iColumn];
1688        double slack = upperSlack_[iColumn]+extra;
1689        deltaSU_[iColumn] = rhsU_[iColumn]-deltaX;
1690        deltaW_[iColumn]=(hHat+wVec_[iColumn]*deltaX)/slack;
1691      }
1692      if (0) {
1693        // different way of calculating
1694        double gamma2 = gamma_*gamma_;
1695        double dZ=0.0;
1696        double dW=0.0;
1697        double zValue = rhsZ_[iColumn];
1698        double gHat = zValue + zVec_[iColumn]*rhsL_[iColumn];
1699        double slackL = lowerSlack_[iColumn]+extra;
1700        double wValue = rhsW_[iColumn];
1701        double hHat = wValue - wVec_[iColumn]*rhsU_[iColumn];
1702        double slackU = upperSlack_[iColumn]+extra;
1703        double q = rhsC_[iColumn]+gamma2 * deltaX +dd;
1704        if (primalR_)
1705          q += deltaX*primalR_[iColumn];
1706        dW = (gHat+hHat -slackL*q + (wValue-zValue)*deltaX)/(slackL+slackU);
1707        dZ = dW + q;
1708        //printf("B %d old %g %g new %g %g\n",iColumn,deltaZ_[iColumn],
1709        //deltaW_[iColumn],dZ,dW);
1710        if (lowerBound(iColumn)) {
1711          if (upperBound(iColumn)) {
1712            //printf("B %d old %g %g new %g %g\n",iColumn,deltaZ_[iColumn],
1713            //deltaW_[iColumn],dZ,dW);
1714            deltaW_[iColumn]=dW;
1715            deltaZ_[iColumn]=dZ;
1716          } else {
1717            // just lower
1718            //printf("L %d old %g new %g\n",iColumn,deltaZ_[iColumn],
1719            //dZ);
1720          }
1721        } else {
1722          assert (upperBound(iColumn));
1723          //printf("U %d old %g new %g\n",iColumn,deltaW_[iColumn],
1724          //dW);
1725        }
1726      }
1727    }
1728  }
1729#if 0
1730  double * check = new double[numberTotal]; 
1731  // Check out rhsC_
1732  multiplyAdd(deltaY_,numberRows_,-1.0,check+numberColumns_,0.0);
1733  CoinZeroN(check,numberColumns_);
1734  matrix_->transposeTimes(1.0,deltaY_,check);
1735  quadraticDjs(check,deltaX_,-1.0);
1736  for (iColumn=0;iColumn<numberTotal;iColumn++) {
1737    check[iColumn] += deltaZ_[iColumn]-deltaW_[iColumn];
1738    if (fabs(check[iColumn]-rhsC_[iColumn])>1.0e-3)
1739      printf("rhsC %d %g %g\n",iColumn,check[iColumn],rhsC_[iColumn]);
1740  }
1741  // Check out rhsZ_
1742  for (iColumn=0;iColumn<numberTotal;iColumn++) {
1743    check[iColumn] += lowerSlack_[iColumn]*deltaZ_[iColumn]+
1744      zVec_[iColumn]*deltaSL_[iColumn];
1745    if (fabs(check[iColumn]-rhsZ_[iColumn])>1.0e-3)
1746      printf("rhsZ %d %g %g\n",iColumn,check[iColumn],rhsZ_[iColumn]);
1747  }
1748  // Check out rhsW_
1749  for (iColumn=0;iColumn<numberTotal;iColumn++) {
1750    check[iColumn] += upperSlack_[iColumn]*deltaW_[iColumn]+
1751      wVec_[iColumn]*deltaSU_[iColumn];
1752    if (fabs(check[iColumn]-rhsW_[iColumn])>1.0e-3)
1753      printf("rhsW %d %g %g\n",iColumn,check[iColumn],rhsW_[iColumn]);
1754  }
1755  delete [] check;
1756#endif
1757  return relativeError;
1758}
1759// createSolution.  Creates solution from scratch
1760int ClpPredictorCorrector::createSolution()
1761{
1762  int numberTotal = numberRows_+numberColumns_;
1763  int iColumn;
1764  double tolerance = primalTolerance();
1765  // See if quadratic objective
1766#ifndef NO_RTTI
1767  ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
1768#else
1769  ClpQuadraticObjective * quadraticObj = NULL;
1770  if (objective_->type()==2)
1771    quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
1772#endif
1773  if (!quadraticObj) {
1774    for (iColumn=0;iColumn<numberTotal;iColumn++) {
1775      if (upper_[iColumn]-lower_[iColumn]>tolerance) 
1776        clearFixed(iColumn);
1777      else
1778        setFixed(iColumn);
1779    }
1780  } else {
1781    // try leaving fixed
1782    for (iColumn=0;iColumn<numberTotal;iColumn++) 
1783      clearFixed(iColumn);
1784  }
1785  double initialValue =1.0e-12;
1786
1787  double maximumObjective=0.0;
1788  double objectiveNorm2=0.0;
1789  getNorms(cost_,numberTotal,maximumObjective,objectiveNorm2);
1790  if (!maximumObjective) {
1791    maximumObjective=1.0; // objective all zero
1792  } 
1793  objectiveNorm2=sqrt(objectiveNorm2)/static_cast<double> (numberTotal);
1794  objectiveNorm_=maximumObjective;
1795  scaleFactor_=1.0;
1796  if (maximumObjective>0.0) {
1797    if (maximumObjective<1.0) {
1798      scaleFactor_=maximumObjective;
1799    } else if (maximumObjective>1.0e4) {
1800      scaleFactor_=maximumObjective/1.0e4;
1801    } 
1802  } 
1803  if (scaleFactor_!=1.0) {
1804    objectiveNorm2*=scaleFactor_;
1805    multiplyAdd(NULL,numberTotal,0.0,cost_,1.0/scaleFactor_);
1806    objectiveNorm_=maximumObjective/scaleFactor_;
1807  }
1808  // See if quadratic objective
1809  if (quadraticObj) {
1810    // If scaled then really scale matrix
1811    double scaleFactor = 
1812      scaleFactor_*optimizationDirection_*objectiveScale_*
1813      rhsScale_;
1814    if ((scalingFlag_>0&&rowScale_)||scaleFactor != 1.0) {
1815      CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
1816      const int * columnQuadratic = quadratic->getIndices();
1817      const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
1818      const int * columnQuadraticLength = quadratic->getVectorLengths();
1819      double * quadraticElement = quadratic->getMutableElements();
1820      int numberColumns = quadratic->getNumCols();
1821      double scale = 1.0/scaleFactor;
1822      if (scalingFlag_>0&&rowScale_) {
1823        for (int iColumn=0;iColumn<numberColumns;iColumn++) {
1824          double scaleI = columnScale_[iColumn]*scale;
1825          for (CoinBigIndex j=columnQuadraticStart[iColumn];
1826               j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
1827            int jColumn = columnQuadratic[j];
1828            double scaleJ = columnScale_[jColumn];
1829            quadraticElement[j] *= scaleI*scaleJ;
1830            objectiveNorm_ = CoinMax(objectiveNorm_,fabs(quadraticElement[j]));
1831          }
1832        }
1833      } else {
1834        // not scaled
1835        for (int iColumn=0;iColumn<numberColumns;iColumn++) {
1836          for (CoinBigIndex j=columnQuadraticStart[iColumn];
1837               j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
1838            quadraticElement[j] *= scale;
1839            objectiveNorm_ = CoinMax(objectiveNorm_,fabs(quadraticElement[j]));
1840          }
1841        }
1842      }
1843    }
1844  }
1845  baseObjectiveNorm_=objectiveNorm_;
1846  //accumulate fixed in dj region (as spare)
1847  //accumulate primal solution in primal region
1848  //DZ in lowerDual
1849  //DW in upperDual
1850  double infiniteCheck=1.0e40;
1851  //double     fakeCheck=1.0e10;
1852  //use deltaX region for work region
1853  for (iColumn=0;iColumn<numberTotal;iColumn++) {
1854    double primalValue = solution_[iColumn];
1855    clearFlagged(iColumn);
1856    clearFixedOrFree(iColumn);
1857    clearLowerBound(iColumn);
1858    clearUpperBound(iColumn);
1859    clearFakeLower(iColumn);
1860    clearFakeUpper(iColumn);
1861    if (!fixed(iColumn)) {
1862      dj_[iColumn]=0.0;
1863      diagonal_[iColumn]=1.0;
1864      deltaX_[iColumn]=1.0;
1865      double lowerValue=lower_[iColumn];
1866      double upperValue=upper_[iColumn];
1867      if (lowerValue>-infiniteCheck) {
1868        if (upperValue<infiniteCheck) {
1869          //upper and lower bounds
1870          setLowerBound(iColumn);
1871          setUpperBound(iColumn);
1872          if (lowerValue>=0.0) {
1873            solution_[iColumn]=lowerValue;
1874          } else if (upperValue<=0.0) {
1875            solution_[iColumn]=upperValue;
1876          } else {
1877            solution_[iColumn]=0.0;
1878          } 
1879        } else {
1880          //just lower bound
1881          setLowerBound(iColumn);
1882          if (lowerValue>=0.0) {
1883            solution_[iColumn]=lowerValue;
1884          } else {
1885            solution_[iColumn]=0.0;
1886          } 
1887        } 
1888      } else {
1889        if (upperValue<infiniteCheck) {
1890          //just upper bound
1891          setUpperBound(iColumn);
1892          if (upperValue<=0.0) {
1893            solution_[iColumn]=upperValue;
1894          } else {
1895            solution_[iColumn]=0.0;
1896          } 
1897        } else {
1898          //free
1899          setFixedOrFree(iColumn);
1900          solution_[iColumn]=0.0;
1901          //std::cout<<" free "<<i<<std::endl;
1902        } 
1903      } 
1904    } else {
1905      setFlagged(iColumn);
1906      setFixedOrFree(iColumn);
1907      setLowerBound(iColumn);
1908      setUpperBound(iColumn);
1909      dj_[iColumn]=primalValue;;
1910      solution_[iColumn]=lower_[iColumn];
1911      diagonal_[iColumn]=0.0;
1912      deltaX_[iColumn]=0.0;
1913    } 
1914  } 
1915  //   modify fixed RHS
1916  multiplyAdd(dj_+numberColumns_,numberRows_,-1.0,rhsFixRegion_,0.0);
1917  //   create plausible RHS?
1918  matrix_->times(-1.0,dj_,rhsFixRegion_);
1919  multiplyAdd(solution_+numberColumns_,numberRows_,1.0,errorRegion_,0.0);
1920  matrix_->times(-1.0,solution_,errorRegion_);
1921  rhsNorm_=maximumAbsElement(errorRegion_,numberRows_);
1922  if (rhsNorm_<1.0) {
1923    rhsNorm_=1.0;
1924  } 
1925  int * rowsDropped = new int [numberRows_];
1926  int returnCode=cholesky_->factorize(diagonal_,rowsDropped);
1927  if (returnCode==-1) {
1928    printf("Out of memory\n");
1929    problemStatus_=4;
1930    return -1;
1931  }
1932  if (cholesky_->status()) {
1933    std::cout << "singular on initial cholesky?" <<std::endl;
1934    cholesky_->resetRowsDropped();
1935    //cholesky_->factorize(rowDropped_);
1936    //if (cholesky_->status()) {
1937      //std::cout << "bad cholesky??? (after retry)" <<std::endl;
1938      //abort();
1939    //}
1940  } 
1941  delete [] rowsDropped;
1942  if (cholesky_->type()<20) {
1943    // not KKT
1944    cholesky_->solve(errorRegion_);
1945    //create information for solution
1946    multiplyAdd(errorRegion_,numberRows_,-1.0,deltaX_+numberColumns_,0.0);
1947    CoinZeroN(deltaX_,numberColumns_);
1948    matrix_->transposeTimes(1.0,errorRegion_,deltaX_);
1949  } else {
1950    // KKT
1951    // reverse sign on solution
1952    multiplyAdd(NULL,numberRows_+numberColumns_,0.0,solution_,-1.0);
1953    solveSystem(deltaX_,errorRegion_,solution_,NULL,NULL,NULL,false);
1954  }
1955  initialValue=1.0e2;
1956  if (rhsNorm_*1.0e-2>initialValue) {
1957    initialValue=rhsNorm_*1.0e-2;
1958  } 
1959  //initialValue = CoinMax(1.0,rhsNorm_);
1960  double smallestBoundDifference=COIN_DBL_MAX;
1961  double * fakeSolution = deltaX_;
1962  for ( iColumn=0;iColumn<numberTotal;iColumn++) {
1963    if (!flagged(iColumn)) {
1964      if (lower_[iColumn]-fakeSolution[iColumn]>initialValue) {
1965        initialValue=lower_[iColumn]-fakeSolution[iColumn];
1966      } 
1967      if (fakeSolution[iColumn]-upper_[iColumn]>initialValue) {
1968        initialValue=fakeSolution[iColumn]-upper_[iColumn];
1969      } 
1970      if (upper_[iColumn]-lower_[iColumn]<smallestBoundDifference) {
1971        smallestBoundDifference=upper_[iColumn]-lower_[iColumn];
1972      } 
1973    } 
1974  } 
1975  solutionNorm_=1.0e-12;
1976  handler_->message(CLP_BARRIER_SAFE,messages_)
1977    <<initialValue<<objectiveNorm_
1978    <<CoinMessageEol;
1979  double extra=1.0e-10;
1980  double largeGap=1.0e15;
1981  //double safeObjectiveValue=2.0*objectiveNorm_;
1982  double safeObjectiveValue=objectiveNorm_+1.0;
1983  double safeFree=1.0e-1*initialValue;
1984  //printf("normal safe dual value of %g, primal value of %g\n",
1985  // safeObjectiveValue,initialValue);
1986  //safeObjectiveValue=CoinMax(2.0,1.0e-1*safeObjectiveValue);
1987  //initialValue=CoinMax(100.0,1.0e-1*initialValue);;
1988  //printf("temp safe dual value of %g, primal value of %g\n",
1989  // safeObjectiveValue,initialValue);
1990  double zwLarge=1.0e2*initialValue;
1991  //zwLarge=1.0e40;
1992  if (cholesky_->choleskyCondition()<0.0&&cholesky_->type()<20) {
1993    // looks bad - play safe
1994    initialValue *=10.0;
1995    safeObjectiveValue *= 10.0;
1996    safeFree *= 10.0;
1997  }
1998  double gamma2 = gamma_*gamma_; // gamma*gamma will be added to diagonal
1999#if 0
2000  fakeSolution[0 ] =   0.072310129 ;
2001  fakeSolution[1 ] =   0.053083871;
2002  fakeSolution[2 ] =      0.178127;
2003  fakeSolution[3 ] =    0.13215151;
2004  fakeSolution[4 ] =   0.072715642;
2005  fakeSolution[5 ] =    0.15680727;
2006  fakeSolution[6 ] =    0.16841689;
2007  fakeSolution[7 ] =   0.093612798 ;
2008  fakeSolution[8 ] =   0.072774891 ;
2009  fakeSolution[9]=1.0;
2010  initialValue=1.0e-5;
2011  safeObjectiveValue=1.0e-5;
2012#endif
2013  // First do primal side
2014  for ( iColumn=0;iColumn<numberTotal;iColumn++) {
2015    if (!flagged(iColumn)) {
2016      double lowerValue=lower_[iColumn];
2017      double upperValue=upper_[iColumn];
2018      double newValue;
2019      double setPrimal=initialValue;
2020      if (quadraticObj) {
2021        // perturb primal solution a bit
2022        //fakeSolution[iColumn]  *= 0.002*CoinDrand48()+0.999;
2023      }
2024      if (lowerBound(iColumn)) {
2025        if (upperBound(iColumn)) {
2026          //upper and lower bounds
2027          if (upperValue-lowerValue>2.0*setPrimal) {
2028            double fakeValue=fakeSolution[iColumn];
2029            if (fakeValue<lowerValue+setPrimal) {
2030              fakeValue=lowerValue+setPrimal;
2031            } 
2032            if (fakeValue>upperValue-setPrimal) {
2033              fakeValue=upperValue-setPrimal;
2034            } 
2035            newValue=fakeValue;
2036          } else {
2037            newValue=0.5*(upperValue+lowerValue);
2038          }
2039        } else {
2040          //just lower bound
2041          double fakeValue=fakeSolution[iColumn];
2042          if (fakeValue<lowerValue+setPrimal) {
2043            fakeValue=lowerValue+setPrimal;
2044          } 
2045          newValue=fakeValue;
2046        } 
2047      } else {
2048        if (upperBound(iColumn)) {
2049          //just upper bound
2050          double fakeValue=fakeSolution[iColumn];
2051          if (fakeValue>upperValue-setPrimal) {
2052            fakeValue=upperValue-setPrimal;
2053          } 
2054          newValue=fakeValue;
2055        } else {
2056          //free
2057          newValue=fakeSolution[iColumn];
2058          if (newValue>=0.0) {
2059            if (newValue<safeFree) {
2060              newValue=safeFree;
2061            } 
2062          } else {
2063            if (newValue>-safeFree) {
2064              newValue=-safeFree;
2065            } 
2066          } 
2067        } 
2068      } 
2069      solution_[iColumn]=newValue;
2070    } else {
2071      // fixed
2072      lowerSlack_[iColumn]=0.0;
2073      upperSlack_[iColumn]=0.0;
2074      solution_[iColumn]=lower_[iColumn];
2075      zVec_[iColumn]=0.0;
2076      wVec_[iColumn]=0.0;
2077      diagonal_[iColumn]=0.0;
2078    } 
2079  }
2080  solutionNorm_ =  maximumAbsElement(solution_,numberTotal);
2081  // Set bounds and do dj including quadratic
2082  largeGap = CoinMax(1.0e7,1.02*solutionNorm_);
2083  CoinPackedMatrix * quadratic = NULL;
2084  const int * columnQuadratic = NULL;
2085  const CoinBigIndex * columnQuadraticStart = NULL;
2086  const int * columnQuadraticLength = NULL;
2087  const double * quadraticElement = NULL;
2088  if (quadraticObj) {
2089    quadratic = quadraticObj->quadraticObjective();
2090    columnQuadratic = quadratic->getIndices();
2091    columnQuadraticStart = quadratic->getVectorStarts();
2092    columnQuadraticLength = quadratic->getVectorLengths();
2093    quadraticElement = quadratic->getElements();
2094  }
2095  double quadraticNorm=0.0;
2096  for ( iColumn=0;iColumn<numberTotal;iColumn++) {
2097    if (!flagged(iColumn)) {
2098      double primalValue=solution_[iColumn];
2099      double lowerValue=lower_[iColumn];
2100      double upperValue=upper_[iColumn];
2101      // Do dj
2102      double reducedCost=cost_[iColumn];
2103      if (lowerBound(iColumn)) {
2104        reducedCost+=linearPerturbation_;
2105      } 
2106      if (upperBound(iColumn)) {
2107        reducedCost-=linearPerturbation_;
2108      }
2109      if (quadraticObj&&iColumn<numberColumns_) {
2110        for (CoinBigIndex j=columnQuadraticStart[iColumn];
2111             j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
2112          int jColumn = columnQuadratic[j];
2113          double valueJ = solution_[jColumn];
2114          double elementValue = quadraticElement[j];
2115          reducedCost += valueJ*elementValue;
2116        }
2117        quadraticNorm = CoinMax(quadraticNorm,fabs(reducedCost));
2118      }
2119      dj_[iColumn]=reducedCost;
2120      if (primalValue>lowerValue+largeGap&&primalValue<upperValue-largeGap) {
2121        clearFixedOrFree(iColumn);
2122        setLowerBound(iColumn);
2123        setUpperBound(iColumn);
2124        lowerValue=CoinMax(lowerValue,primalValue-largeGap);
2125        upperValue=CoinMin(upperValue,primalValue+largeGap);
2126        lower_[iColumn]=lowerValue;
2127        upper_[iColumn]=upperValue;
2128      }
2129    }
2130  }
2131  safeObjectiveValue=CoinMax(safeObjectiveValue,quadraticNorm);
2132  for ( iColumn=0;iColumn<numberTotal;iColumn++) {
2133    if (!flagged(iColumn)) {
2134      double primalValue=solution_[iColumn];
2135      double lowerValue=lower_[iColumn];
2136      double upperValue=upper_[iColumn];
2137      double reducedCost=dj_[iColumn];
2138      double low=0.0;
2139      double high=0.0;
2140      if (lowerBound(iColumn)) {
2141        if (upperBound(iColumn)) {
2142          //upper and lower bounds
2143          if (upperValue-lowerValue>2.0*initialValue) {
2144            low=primalValue-lowerValue;
2145            high=upperValue-primalValue;
2146          } else {
2147            low=initialValue;
2148            high=initialValue;
2149          }
2150          double s = low+extra;
2151          double ratioZ;
2152          if (s<zwLarge) {
2153            ratioZ=1.0;
2154          } else {
2155            ratioZ=sqrt(zwLarge/s);
2156          } 
2157          double t = high+extra;
2158          double ratioT;
2159          if (t<zwLarge) {
2160            ratioT=1.0;
2161          } else {
2162            ratioT=sqrt(zwLarge/t);
2163          } 
2164          //modify s and t
2165          if (s>largeGap) {
2166            s=largeGap;
2167          } 
2168          if (t>largeGap) {
2169            t=largeGap;
2170          } 
2171          //modify if long long way away from bound
2172          if (reducedCost>=0.0) {
2173            zVec_[iColumn]=reducedCost + safeObjectiveValue*ratioZ;
2174            zVec_[iColumn]=CoinMax(reducedCost, safeObjectiveValue*ratioZ);
2175            wVec_[iColumn]=safeObjectiveValue*ratioT;
2176          } else {
2177            zVec_[iColumn]=safeObjectiveValue*ratioZ;
2178            wVec_[iColumn]=-reducedCost + safeObjectiveValue*ratioT;
2179            wVec_[iColumn]=CoinMax(-reducedCost , safeObjectiveValue*ratioT);
2180          }
2181          double gammaTerm = gamma2;
2182          if (primalR_)
2183            gammaTerm += primalR_[iColumn];
2184          diagonal_[iColumn] = (t*s)/
2185            (s*wVec_[iColumn]+t*zVec_[iColumn]+gammaTerm*t*s);
2186        } else {
2187          //just lower bound
2188          low=primalValue-lowerValue;
2189          high=0.0;
2190          double s = low+extra;
2191          double ratioZ;
2192          if (s<zwLarge) {
2193            ratioZ=1.0;
2194          } else {
2195            ratioZ=sqrt(zwLarge/s);
2196          } 
2197          //modify s
2198          if (s>largeGap) {
2199            s=largeGap;
2200          } 
2201          if (reducedCost>=0.0) {
2202            zVec_[iColumn]=reducedCost + safeObjectiveValue*ratioZ;
2203            zVec_[iColumn]=CoinMax(reducedCost , safeObjectiveValue*ratioZ);
2204            wVec_[iColumn]=0.0;
2205          } else {
2206            zVec_[iColumn]=safeObjectiveValue*ratioZ;
2207            wVec_[iColumn]=0.0;
2208          } 
2209          double gammaTerm = gamma2;
2210          if (primalR_)
2211            gammaTerm += primalR_[iColumn];
2212          diagonal_[iColumn] = s/(zVec_[iColumn]+s*gammaTerm);
2213        } 
2214      } else {
2215        if (upperBound(iColumn)) {
2216          //just upper bound
2217          low=0.0;
2218          high=upperValue-primalValue;
2219          double t = high+extra;
2220          double ratioT;
2221          if (t<zwLarge) {
2222            ratioT=1.0;
2223          } else {
2224            ratioT=sqrt(zwLarge/t);
2225          } 
2226          //modify t
2227          if (t>largeGap) {
2228            t=largeGap;
2229          } 
2230          if (reducedCost>=0.0) {
2231            zVec_[iColumn]=0.0;
2232            wVec_[iColumn]=safeObjectiveValue*ratioT;
2233          } else {
2234            zVec_[iColumn]=0.0;
2235            wVec_[iColumn]=-reducedCost + safeObjectiveValue*ratioT;
2236            wVec_[iColumn]=CoinMax(-reducedCost , safeObjectiveValue*ratioT);
2237          } 
2238          double gammaTerm = gamma2;
2239          if (primalR_)
2240            gammaTerm += primalR_[iColumn];
2241          diagonal_[iColumn] =  t/(wVec_[iColumn]+t*gammaTerm);
2242        } 
2243      } 
2244      lowerSlack_[iColumn]=low;
2245      upperSlack_[iColumn]=high;
2246    }
2247  }
2248#if 0
2249  if (solution_[0]>0.0) {
2250    for (int i=0;i<numberTotal;i++)
2251      printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n",i,fabs(solution_[i]),
2252             diagonal_[i],fabs(dj_[i]),
2253             lowerSlack_[i],zVec_[i],
2254             upperSlack_[i],wVec_[i]);
2255  } else {
2256    for (int i=0;i<numberTotal;i++)
2257      printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n",i,fabs(solution_[i]),
2258             diagonal_[i],fabs(dj_[i]),
2259             upperSlack_[i],wVec_[i],
2260             lowerSlack_[i],zVec_[i] );
2261  }
2262  exit(66);
2263#endif
2264  return 0;
2265}
2266// complementarityGap.  Computes gap
2267//phase 0=as is , 1 = after predictor , 2 after corrector
2268double ClpPredictorCorrector::complementarityGap(int & numberComplementarityPairs,
2269                                                 int & numberComplementarityItems,
2270                                                 const int phase)
2271{
2272  double gap=0.0;
2273  //seems to be same coding for phase = 1 or 2
2274  numberComplementarityPairs=0;
2275  numberComplementarityItems=0;
2276  int numberTotal = numberRows_+numberColumns_;
2277  double toleranceGap=0.0;
2278  double largestGap=0.0;
2279  double smallestGap=COIN_DBL_MAX;
2280  //seems to be same coding for phase = 1 or 2
2281  int numberNegativeGaps=0;
2282  double sumNegativeGap=0.0;
2283  double largeGap=1.0e2*solutionNorm_;
2284  if (largeGap<1.0e10) {
2285    largeGap=1.0e10;
2286  } 
2287  largeGap=1.0e30;
2288  double dualTolerance =  dblParam_[ClpDualTolerance];
2289  double primalTolerance =  dblParam_[ClpPrimalTolerance];
2290  dualTolerance=dualTolerance/scaleFactor_;
2291  for (int iColumn=0;iColumn<numberTotal;iColumn++) {
2292    if (!fixedOrFree(iColumn)) {
2293      numberComplementarityPairs++;
2294      //can collapse as if no lower bound both zVec and deltaZ 0.0
2295      double newZ=0.0;
2296      double newW=0.0;
2297      if (lowerBound(iColumn)) {
2298        numberComplementarityItems++;
2299        double dualValue;
2300        double primalValue;
2301        if (!phase) {
2302          dualValue=zVec_[iColumn];
2303          primalValue=lowerSlack_[iColumn];
2304        } else {
2305          double change;
2306          change =solution_[iColumn]+deltaX_[iColumn]-lowerSlack_[iColumn]-lower_[iColumn];
2307          dualValue=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn];
2308          newZ=dualValue;
2309          primalValue=lowerSlack_[iColumn]+actualPrimalStep_*change;
2310        }
2311        //reduce primalValue
2312        if (primalValue>largeGap) {
2313          primalValue=largeGap;
2314        } 
2315        double gapProduct=dualValue*primalValue;
2316        if (gapProduct<0.0) {
2317          //cout<<"negative gap component "<<iColumn<<" "<<dualValue<<" "<<
2318              //primalValue<<endl;
2319          numberNegativeGaps++;
2320          sumNegativeGap-=gapProduct;
2321          gapProduct=0.0;
2322        } 
2323        gap+=gapProduct;
2324        //printf("l %d prim %g dual %g totalGap %g\n",
2325        //   iColumn,primalValue,dualValue,gap);
2326        if (gapProduct>largestGap) {
2327          largestGap=gapProduct;
2328        }
2329        smallestGap = CoinMin(smallestGap,gapProduct);
2330        if (dualValue>dualTolerance&&primalValue>primalTolerance) {
2331          toleranceGap+=dualValue*primalValue;
2332        } 
2333      } 
2334      if (upperBound(iColumn)) {
2335        numberComplementarityItems++;
2336        double dualValue;
2337        double primalValue;
2338        if (!phase) {
2339          dualValue=wVec_[iColumn];
2340          primalValue=upperSlack_[iColumn];
2341        } else {
2342          double change;
2343          change =upper_[iColumn]-solution_[iColumn]-deltaX_[iColumn]-upperSlack_[iColumn];
2344          dualValue=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn];
2345          newW=dualValue;
2346          primalValue=upperSlack_[iColumn]+actualPrimalStep_*change;
2347        } 
2348        //reduce primalValue
2349        if (primalValue>largeGap) {
2350          primalValue=largeGap;
2351        } 
2352        double gapProduct=dualValue*primalValue;
2353        if (gapProduct<0.0) {
2354          //cout<<"negative gap component "<<iColumn<<" "<<dualValue<<" "<<
2355              //primalValue<<endl;
2356          numberNegativeGaps++;
2357          sumNegativeGap-=gapProduct;
2358          gapProduct=0.0;
2359        } 
2360        gap+=gapProduct;
2361        //printf("u %d prim %g dual %g totalGap %g\n",
2362        //   iColumn,primalValue,dualValue,gap);
2363        if (gapProduct>largestGap) {
2364          largestGap=gapProduct;
2365        } 
2366        if (dualValue>dualTolerance&&primalValue>primalTolerance) {
2367          toleranceGap+=dualValue*primalValue;
2368        } 
2369      } 
2370    } 
2371  }
2372  //if (numberIterations_>4)
2373  //exit(9);
2374  if (!phase&&numberNegativeGaps) {
2375      handler_->message(CLP_BARRIER_NEGATIVE_GAPS,messages_)
2376    <<numberNegativeGaps<<sumNegativeGap
2377    <<CoinMessageEol;
2378  } 
2379 
2380  //in case all free!
2381  if (!numberComplementarityPairs) {
2382    numberComplementarityPairs=1;
2383  } 
2384#ifdef SOME_DEBUG
2385  printf("with d,p steps %g,%g gap %g - smallest %g, largest %g, pairs %d\n",
2386         actualDualStep_,actualPrimalStep_,
2387         gap,smallestGap,largestGap,numberComplementarityPairs);
2388#endif
2389  return gap;
2390}
2391// setupForSolve.
2392//phase 0=affine , 1 = corrector , 2 = primal-dual
2393void ClpPredictorCorrector::setupForSolve(const int phase)
2394{
2395  double extra =eExtra;
2396  int numberTotal = numberRows_ + numberColumns_;
2397  int iColumn;
2398#ifdef SOME_DEBUG
2399  printf("phase %d in setupForSolve, mu %.18g\n",phase,mu_);
2400#endif
2401  double gamma2 = gamma_*gamma_; // gamma*gamma will be added to diagonal
2402  switch (phase) {
2403  case 0:
2404    CoinMemcpyN(errorRegion_,numberRows_,rhsB_);
2405    if (delta_||dualR_) {
2406      // add in regularization
2407      double delta2 = delta_*delta_;
2408      for (int iRow=0;iRow<numberRows_;iRow++) {
2409        rhsB_[iRow] -= delta2*dual_[iRow];
2410        if (dualR_)
2411          rhsB_[iRow] -= dualR_[iRow]*dual_[iRow];
2412      }
2413    }
2414    for (iColumn=0;iColumn<numberTotal;iColumn++) {
2415      rhsC_[iColumn]=0.0;
2416      rhsU_[iColumn]=0.0;
2417      rhsL_[iColumn]=0.0;
2418      rhsZ_[iColumn]=0.0;
2419      rhsW_[iColumn]=0.0;
2420      if (!flagged(iColumn)) {
2421        rhsC_[iColumn] = dj_[iColumn]-zVec_[iColumn]+wVec_[iColumn];
2422        rhsC_[iColumn] += gamma2*solution_[iColumn];
2423        if (primalR_)
2424          rhsC_[iColumn] += primalR_[iColumn]*solution_[iColumn];
2425        if (lowerBound(iColumn)) {
2426          rhsZ_[iColumn] = -zVec_[iColumn]*(lowerSlack_[iColumn]+extra);
2427          rhsL_[iColumn] = CoinMax(0.0,(lower_[iColumn]+lowerSlack_[iColumn])-solution_[iColumn]);
2428        } 
2429        if (upperBound(iColumn)) {
2430          rhsW_[iColumn] = -wVec_[iColumn]*(upperSlack_[iColumn]+extra);
2431          rhsU_[iColumn] = CoinMin(0.0,(upper_[iColumn]-upperSlack_[iColumn])-solution_[iColumn]);
2432        }
2433      }
2434    } 
2435    if (0) {
2436      int i=1324;
2437      printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n",i,solution_[i],
2438             diagonal_[i],dj_[i],
2439             lowerSlack_[i],zVec_[i],
2440             upperSlack_[i],wVec_[i]);
2441      printf("%d %.18g %.18g %.18g %.18g %.18g\n",i,rhsC_[i],
2442             rhsZ_[i],rhsL_[i],
2443             rhsW_[i],rhsU_[i]);
2444    }
2445#if 0
2446    for (int i=0;i<3;i++) {
2447      if (!fabs(rhsZ_[i]))
2448        rhsZ_[i]=0.0;
2449      if (!fabs(rhsW_[i]))
2450        rhsW_[i]=0.0;
2451      if (!fabs(rhsU_[i]))
2452        rhsU_[i]=0.0;
2453      if (!fabs(rhsL_[i]))
2454        rhsL_[i]=0.0;
2455    }
2456    if (solution_[0]>0.0) {
2457      for (int i=0;i<3;i++)
2458        printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n",i,solution_[i],
2459               diagonal_[i],dj_[i],
2460               lowerSlack_[i],zVec_[i],
2461               upperSlack_[i],wVec_[i]);
2462      for (int i=0;i<3;i++)
2463        printf("%d %.18g %.18g %.18g %.18g %.18g\n",i,rhsC_[i],
2464               rhsZ_[i],rhsL_[i],
2465               rhsW_[i],rhsU_[i]);
2466    } else {
2467      for (int i=0;i<3;i++)
2468        printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n",i,solution_[i],
2469               diagonal_[i],dj_[i],
2470               lowerSlack_[i],zVec_[i],
2471               upperSlack_[i],wVec_[i]);
2472      for (int i=0;i<3;i++)
2473        printf("%d %.18g %.18g %.18g %.18g %.18g\n",i,rhsC_[i],
2474               rhsZ_[i],rhsL_[i],
2475               rhsW_[i],rhsU_[i]);
2476    }
2477#endif
2478    break;
2479  case 1:
2480    // could be stored in delta2?
2481    for (iColumn=0;iColumn<numberTotal;iColumn++) {
2482      rhsZ_[iColumn]=0.0;
2483      rhsW_[iColumn]=0.0;
2484      if (!flagged(iColumn)) {
2485        if (lowerBound(iColumn)) {
2486          rhsZ_[iColumn] = mu_ -zVec_[iColumn]*(lowerSlack_[iColumn]+extra)
2487            - deltaZ_[iColumn]*deltaX_[iColumn];
2488          // To bring in line with OSL
2489          rhsZ_[iColumn] += deltaZ_[iColumn]*rhsL_[iColumn];
2490        } 
2491        if (upperBound(iColumn)) {
2492          rhsW_[iColumn] = mu_ -wVec_[iColumn]*(upperSlack_[iColumn]+extra)
2493            +deltaW_[iColumn]*deltaX_[iColumn];
2494          // To bring in line with OSL
2495          rhsW_[iColumn] -= deltaW_[iColumn]*rhsU_[iColumn];
2496        } 
2497      } 
2498    } 
2499#if 0
2500    for (int i=0;i<numberTotal;i++) {
2501      if (!fabs(rhsZ_[i]))
2502        rhsZ_[i]=0.0;
2503      if (!fabs(rhsW_[i]))
2504        rhsW_[i]=0.0;
2505      if (!fabs(rhsU_[i]))
2506        rhsU_[i]=0.0;
2507      if (!fabs(rhsL_[i]))
2508        rhsL_[i]=0.0;
2509    }
2510    if (solution_[0]>0.0) {
2511      for (int i=0;i<numberTotal;i++)
2512        printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n",i,fabs(solution_[i]),
2513               diagonal_[i],fabs(dj_[i]),
2514               lowerSlack_[i],zVec_[i],
2515               upperSlack_[i],wVec_[i]);
2516      for (int i=0;i<numberTotal;i++)
2517        printf("%d %.18g %.18g %.18g %.18g %.18g\n",i,fabs(rhsC_[i]),
2518               rhsZ_[i],rhsL_[i],
2519               rhsW_[i],rhsU_[i]);
2520    } else {
2521      for (int i=0;i<numberTotal;i++)
2522        printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n",i,fabs(solution_[i]),
2523               diagonal_[i],fabs(dj_[i]),
2524               upperSlack_[i],wVec_[i],
2525               lowerSlack_[i],zVec_[i] );
2526      for (int i=0;i<numberTotal;i++)
2527        printf("%d %.18g %.18g %.18g %.18g %.18g\n",i,fabs(rhsC_[i]),
2528               rhsW_[i],rhsU_[i],
2529               rhsZ_[i],rhsL_[i]);
2530    }
2531    exit(66);
2532#endif
2533    break;
2534  case 2:
2535    CoinMemcpyN(errorRegion_,numberRows_,rhsB_);
2536    for (iColumn=0;iColumn<numberTotal;iColumn++) {
2537      rhsZ_[iColumn]=0.0;
2538      rhsW_[iColumn]=0.0;
2539      if (!flagged(iColumn)) {
2540        if (lowerBound(iColumn)) {
2541          rhsZ_[iColumn] = mu_ - zVec_[iColumn]*(lowerSlack_[iColumn]+extra);
2542        } 
2543        if (upperBound(iColumn)) {
2544          rhsW_[iColumn] = mu_ - wVec_[iColumn]*(upperSlack_[iColumn]+extra);
2545        }
2546      }
2547    } 
2548    break;
2549  case 3:
2550    {
2551      double minBeta = 0.1*mu_;
2552      double maxBeta = 10.0*mu_;
2553      double dualStep = CoinMin(1.0,actualDualStep_+0.1);
2554      double primalStep = CoinMin(1.0,actualPrimalStep_+0.1);
2555#ifdef SOME_DEBUG
2556      printf("good complementarity range %g to %g\n",minBeta,maxBeta);
2557#endif
2558      //minBeta=0.0;
2559      //maxBeta=COIN_DBL_MAX;
2560      for (iColumn=0;iColumn<numberTotal;iColumn++) {
2561        if (!flagged(iColumn)) {
2562          if (lowerBound(iColumn)) {
2563            double change = -rhsL_[iColumn] + deltaX_[iColumn];
2564            double dualValue=zVec_[iColumn]+dualStep*deltaZ_[iColumn];
2565            double primalValue=lowerSlack_[iColumn]+primalStep*change;
2566            double gapProduct=dualValue*primalValue;
2567            if (gapProduct>0.0&&dualValue<0.0)
2568              gapProduct = - gapProduct;
2569#ifdef FULL_DEBUG
2570            delta2Z_[iColumn]=gapProduct;
2571            if (delta2Z_[iColumn]<minBeta||delta2Z_[iColumn]>maxBeta)
2572              printf("lower %d primal %g, dual %g, gap %g\n",
2573                     iColumn,primalValue,dualValue,gapProduct);
2574#endif
2575            double value= 0.0;
2576            if (gapProduct<minBeta) {
2577              value= 2.0*(minBeta-gapProduct);
2578              value= (mu_-gapProduct);
2579              value= (minBeta-gapProduct);
2580              assert (value>0.0);
2581            } else if (gapProduct>maxBeta) {
2582              value= CoinMax(maxBeta-gapProduct,-maxBeta);
2583              assert (value<0.0);
2584            }
2585            rhsZ_[iColumn] += value;
2586          } 
2587          if (upperBound(iColumn)) {
2588            double change = rhsU_[iColumn]-deltaX_[iColumn];
2589            double dualValue=wVec_[iColumn]+dualStep*deltaW_[iColumn];
2590            double primalValue=upperSlack_[iColumn]+primalStep*change;
2591            double gapProduct=dualValue*primalValue;
2592            if (gapProduct>0.0&&dualValue<0.0)
2593              gapProduct = - gapProduct;
2594#ifdef FULL_DEBUG
2595            delta2W_[iColumn]=gapProduct;
2596            if (delta2W_[iColumn]<minBeta||delta2W_[iColumn]>maxBeta)
2597              printf("upper %d primal %g, dual %g, gap %g\n",
2598                     iColumn,primalValue,dualValue,gapProduct);
2599#endif
2600            double value= 0.0;
2601            if (gapProduct<minBeta) {
2602              value= (minBeta-gapProduct);
2603              assert (value>0.0);
2604            } else if (gapProduct>maxBeta) {
2605              value= CoinMax(maxBeta-gapProduct,-maxBeta);
2606              assert (value<0.0);
2607            }
2608            rhsW_[iColumn] += value;
2609          } 
2610        } 
2611      }
2612    }
2613    break;
2614  } /* endswitch */
2615  if (cholesky_->type()<20) {
2616    for (iColumn=0;iColumn<numberTotal;iColumn++) {
2617      double value = rhsC_[iColumn];
2618      double zValue = rhsZ_[iColumn];
2619      double wValue = rhsW_[iColumn];
2620#if 0
2621#if 1
2622      if (phase==0) {
2623        // more accurate
2624        value = dj[iColumn];
2625        zValue=0.0;
2626        wValue=0.0;
2627      } else if (phase==2) {
2628        // more accurate
2629        value = dj[iColumn];
2630        zValue=mu_;
2631        wValue=mu_;
2632      }
2633#endif
2634      assert (rhsL_[iColumn]>=0.0);
2635      assert (rhsU_[iColumn]<=0.0);
2636      if (lowerBound(iColumn)) {
2637        value += (-zVec_[iColumn]*rhsL_[iColumn]-zValue)/
2638          (lowerSlack_[iColumn]+extra);
2639      }
2640      if (upperBound(iColumn)) {
2641        value += (wValue-wVec_[iColumn]*rhsU_[iColumn])/
2642          (upperSlack_[iColumn]+extra);
2643      }
2644#else
2645      if (lowerBound(iColumn)) {
2646        double gHat = zValue + zVec_[iColumn]*rhsL_[iColumn];
2647        value -= gHat/(lowerSlack_[iColumn]+extra);
2648      }
2649      if (upperBound(iColumn)) {
2650        double hHat = wValue - wVec_[iColumn]*rhsU_[iColumn];
2651        value += hHat/(upperSlack_[iColumn]+extra);
2652      }
2653#endif
2654      workArray_[iColumn]=diagonal_[iColumn]*value;
2655    } 
2656#if 0
2657    if (solution_[0]>0.0) {
2658      for (int i=0;i<numberTotal;i++)
2659        printf("%d %.18g\n",i,workArray_[i]);
2660    } else {
2661      for (int i=0;i<numberTotal;i++)
2662        printf("%d %.18g\n",i,workArray_[i]);
2663    }
2664    exit(66);
2665#endif
2666  } else {
2667    // KKT
2668    for (iColumn=0;iColumn<numberTotal;iColumn++) {
2669      double value = rhsC_[iColumn];
2670      double zValue = rhsZ_[iColumn];
2671      double wValue = rhsW_[iColumn];
2672      if (lowerBound(iColumn)) {
2673        double gHat = zValue + zVec_[iColumn]*rhsL_[iColumn];
2674        value -= gHat/(lowerSlack_[iColumn]+extra);
2675      }
2676      if (upperBound(iColumn)) {
2677        double hHat = wValue - wVec_[iColumn]*rhsU_[iColumn];
2678        value += hHat/(upperSlack_[iColumn]+extra);
2679      }
2680      workArray_[iColumn]=value;
2681    }
2682  }
2683}
2684//method: sees if looks plausible change in complementarity
2685bool ClpPredictorCorrector::checkGoodMove(const bool doCorrector,
2686                                          double & bestNextGap,
2687                                          bool allowIncreasingGap)
2688{
2689  const double beta3 = 0.99997;
2690  bool goodMove=false;
2691  int nextNumber;
2692  int nextNumberItems;
2693  int numberTotal = numberRows_+numberColumns_;
2694  double returnGap=bestNextGap;
2695  double nextGap=complementarityGap(nextNumber,nextNumberItems,2);
2696#ifndef NO_RTTI
2697  ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
2698#else
2699  ClpQuadraticObjective * quadraticObj = NULL;
2700  if (objective_->type()==2)
2701    quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
2702#endif
2703  if (nextGap>bestNextGap&&nextGap>0.9*complementarityGap_&&doCorrector
2704      &&!quadraticObj&&!allowIncreasingGap) {
2705#ifdef SOME_DEBUG
2706    printf("checkGood phase 1 next gap %.18g, phase 0 %.18g, old gap %.18g\n",
2707           nextGap,bestNextGap,complementarityGap_);
2708#endif
2709    return false;
2710  } else {
2711    returnGap=nextGap;
2712  }
2713  double step;
2714  if (actualDualStep_>actualPrimalStep_) {
2715    step=actualDualStep_;
2716  } else {
2717    step=actualPrimalStep_;
2718  } 
2719  double testValue=1.0-step*(1.0-beta3);
2720  //testValue=0.0;
2721  testValue*=complementarityGap_;
2722  if (nextGap<testValue) {
2723    //std::cout <<"predicted duality gap "<<nextGap<<std::endl;
2724    goodMove=true;
2725  } else if(doCorrector) {
2726    //if (actualDualStep_<actualPrimalStep_) {
2727      //step=actualDualStep_;
2728    //} else {
2729      //step=actualPrimalStep_;
2730    //}
2731    double gap = bestNextGap;
2732    goodMove=checkGoodMove2(step,gap,allowIncreasingGap);
2733    if (goodMove)
2734      returnGap=gap;
2735  } else {
2736    goodMove=true;
2737  } 
2738  if (goodMove)
2739    goodMove=checkGoodMove2(step,bestNextGap,allowIncreasingGap);
2740  // Say good if small
2741  //if (quadraticObj) {
2742  if (CoinMax(actualDualStep_,actualPrimalStep_)<1.0e-6)
2743    goodMove=true;
2744  if (!goodMove) {
2745    //try smaller of two
2746    if (actualDualStep_<actualPrimalStep_) {
2747      step=actualDualStep_;
2748    } else {
2749      step=actualPrimalStep_;
2750    } 
2751    if (step>1.0) {
2752      step=1.0;
2753    } 
2754    actualPrimalStep_=step;
2755    //if (quadraticObj)
2756    //actualPrimalStep_ *=0.5;
2757    actualDualStep_=step;
2758    goodMove=checkGoodMove2(step,bestNextGap,allowIncreasingGap);
2759    int pass=0;
2760    while (!goodMove) {
2761      pass++;
2762      double gap = bestNextGap;
2763      goodMove=checkGoodMove2(step,gap,allowIncreasingGap);
2764      if (goodMove||pass>3) {
2765        returnGap=gap;
2766        break;
2767      }
2768      if (step<1.0e-4) {
2769        break;
2770      } 
2771      step*=0.5;
2772      actualPrimalStep_=step;
2773      //if (quadraticObj)
2774      //actualPrimalStep_ *=0.5;
2775      actualDualStep_=step;
2776    } /* endwhile */
2777    if (doCorrector) {
2778      //say bad move if both small
2779      if (numberIterations_&1) {
2780        if (actualPrimalStep_<1.0e-2&&actualDualStep_<1.0e-2) {
2781          goodMove=false;
2782        } 
2783      } else {
2784        if (actualPrimalStep_<1.0e-5&&actualDualStep_<1.0e-5) {
2785          goodMove=false;
2786        } 
2787        if (actualPrimalStep_*actualDualStep_<1.0e-20) {
2788          goodMove=false;
2789        } 
2790      } 
2791    } 
2792  } 
2793  if (goodMove) {
2794    //compute delta in objectives
2795    double deltaObjectivePrimal=0.0;
2796    double deltaObjectiveDual=
2797      innerProduct(deltaY_,numberRows_,
2798                   rhsFixRegion_);
2799    double error=0.0;
2800    double * workArray = workArray_;
2801    CoinZeroN(workArray,numberColumns_);
2802    CoinMemcpyN(deltaY_,numberRows_,workArray+numberColumns_);
2803    matrix_->transposeTimes(-1.0,deltaY_,workArray);
2804    //double sumPerturbCost=0.0;
2805    for (int iColumn=0;iColumn<numberTotal;iColumn++) {
2806      if (!flagged(iColumn)) {
2807        if (lowerBound(iColumn)) {
2808          //sumPerturbCost+=deltaX_[iColumn];
2809          deltaObjectiveDual+=deltaZ_[iColumn]*lower_[iColumn];
2810        } 
2811        if (upperBound(iColumn)) {
2812          //sumPerturbCost-=deltaX_[iColumn];
2813          deltaObjectiveDual-=deltaW_[iColumn]*upper_[iColumn];
2814        } 
2815        double change = fabs(workArray_[iColumn]-deltaZ_[iColumn]+deltaW_[iColumn]);
2816        error = CoinMax(change,error);
2817      } 
2818      deltaObjectivePrimal += cost_[iColumn] * deltaX_[iColumn];
2819    } 
2820    //deltaObjectivePrimal+=sumPerturbCost*linearPerturbation_;
2821    double testValue;
2822    if (error>0.0) {
2823      testValue=1.0e1*CoinMax(maximumDualError_,1.0e-12)/error;
2824    } else {
2825      testValue=1.0e1;
2826    } 
2827    // If quadratic then primal step may compensate
2828    if (testValue<actualDualStep_&&!quadraticObj) {
2829      handler_->message(CLP_BARRIER_REDUCING,messages_)
2830      <<"dual"<<actualDualStep_
2831      << testValue
2832      <<CoinMessageEol;
2833      actualDualStep_=testValue;
2834    } 
2835  } 
2836  if (maximumRHSError_<1.0e1*solutionNorm_*primalTolerance()
2837                            &&maximumRHSChange_>1.0e-16*solutionNorm_) {
2838    //check change in AX not too much
2839    //??? could be dropped row going infeasible
2840    double ratio = 1.0e1*CoinMax(maximumRHSError_,1.0e-12)/maximumRHSChange_;
2841    if (ratio<actualPrimalStep_) {
2842      handler_->message(CLP_BARRIER_REDUCING,messages_)
2843      <<"primal"<<actualPrimalStep_
2844      <<ratio
2845      <<CoinMessageEol;
2846      if (ratio>1.0e-6) {
2847        actualPrimalStep_=ratio;
2848      } else {
2849        actualPrimalStep_=ratio;
2850        //std::cout <<"sign we should be stopping"<<std::endl;
2851      } 
2852    } 
2853  }
2854  if (goodMove)
2855    bestNextGap=returnGap;
2856  return goodMove;
2857}
2858//:  checks for one step size
2859bool ClpPredictorCorrector::checkGoodMove2(double move,
2860                                           double & bestNextGap,
2861                                           bool allowIncreasingGap)
2862{
2863  double complementarityMultiplier =1.0/numberComplementarityPairs_;
2864  const double gamma = 1.0e-8;
2865  const double gammap = 1.0e-8;
2866  double gammad = 1.0e-8;
2867  int nextNumber;
2868  int nextNumberItems;
2869  double nextGap=complementarityGap(nextNumber,nextNumberItems,2);
2870  if (nextGap>bestNextGap&&!allowIncreasingGap)
2871    return false;
2872  double lowerBoundGap = gamma*nextGap*complementarityMultiplier;
2873  bool goodMove=true;
2874  int iColumn;
2875  for ( iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
2876    if (!flagged(iColumn)) {
2877      if (lowerBound(iColumn)) {
2878        double part1=lowerSlack_[iColumn]+actualPrimalStep_*deltaSL_[iColumn];
2879        double part2=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn];
2880        if (part1*part2<lowerBoundGap) {
2881          goodMove=false;
2882          break;
2883        } 
2884      } 
2885      if (upperBound(iColumn)) {
2886        double part1=upperSlack_[iColumn]+actualPrimalStep_*deltaSU_[iColumn];
2887        double part2=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn];
2888        if (part1*part2<lowerBoundGap) {
2889          goodMove=false;
2890          break;
2891        } 
2892      } 
2893    } 
2894  } 
2895   double * nextDj=NULL;
2896   double maximumDualError = maximumDualError_;
2897#ifndef NO_RTTI
2898  ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
2899#else
2900  ClpQuadraticObjective * quadraticObj = NULL;
2901  if (objective_->type()==2)
2902    quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
2903#endif
2904  if (quadraticObj) {
2905    // change gammad
2906    gammad=1.0e-4;
2907    double gamma2 = gamma_*gamma_;
2908    nextDj = new double [numberColumns_];
2909    double * nextSolution = new double [numberColumns_];
2910    // put next primal into nextSolution
2911    for ( iColumn=0;iColumn<numberColumns_;iColumn++) {
2912      if (!flagged(iColumn)) {
2913        nextSolution[iColumn]=solution_[iColumn]+
2914          actualPrimalStep_*deltaX_[iColumn];
2915      } else {
2916        nextSolution[iColumn]=solution_[iColumn];
2917      }
2918    }
2919    // do reduced costs
2920    CoinMemcpyN(cost_,numberColumns_,nextDj);
2921    matrix_->transposeTimes(-1.0,dual_,nextDj);
2922    matrix_->transposeTimes(-actualDualStep_,deltaY_,nextDj);
2923    quadraticDjs(nextDj,nextSolution,1.0);
2924    delete [] nextSolution;
2925    CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
2926    const int * columnQuadraticLength = quadratic->getVectorLengths();
2927    for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
2928      if (!fixedOrFree(iColumn)) {
2929        double newZ=0.0;
2930        double newW=0.0;
2931        if (lowerBound(iColumn)) {
2932          newZ=zVec_[iColumn]+actualDualStep_*deltaZ_[iColumn];
2933        } 
2934        if (upperBound(iColumn)) {
2935          newW=wVec_[iColumn]+actualDualStep_*deltaW_[iColumn];
2936        } 
2937        if (columnQuadraticLength[iColumn]) {
2938          double gammaTerm = gamma2;
2939          if (primalR_)
2940            gammaTerm += primalR_[iColumn];
2941          //double dualInfeasibility=
2942          //dj_[iColumn]-zVec_[iColumn]+wVec_[iColumn]
2943          //+gammaTerm*solution_[iColumn];
2944          double newInfeasibility=
2945            nextDj[iColumn]-newZ+newW
2946            +gammaTerm*(solution_[iColumn]+actualPrimalStep_*deltaX_[iColumn]);
2947          maximumDualError = CoinMax(maximumDualError,newInfeasibility);
2948          //if (fabs(newInfeasibility)>CoinMax(2000.0*maximumDualError_,1.0e-2)) {
2949          //if (dualInfeasibility*newInfeasibility<0.0) {
2950          //  printf("%d current %g next %g\n",iColumn,dualInfeasibility,
2951          //       newInfeasibility);
2952          //  goodMove=false;
2953          //}
2954          //}
2955        }
2956      } 
2957    }
2958    delete [] nextDj;
2959  }
2960 //      Satisfy g_p(alpha)?
2961  if (rhsNorm_>solutionNorm_) {
2962    solutionNorm_=rhsNorm_;
2963  } 
2964  double errorCheck=maximumRHSError_/solutionNorm_;
2965  if (errorCheck<maximumBoundInfeasibility_) {
2966    errorCheck=maximumBoundInfeasibility_;
2967  } 
2968  // scale back move
2969  move = CoinMin(move,0.95);
2970  //scale
2971  if ((1.0-move)*errorCheck>primalTolerance()) {
2972    if (nextGap<gammap*(1.0-move)*errorCheck) {
2973      goodMove=false;
2974    } 
2975  } 
2976  //      Satisfy g_d(alpha)?
2977  errorCheck=maximumDualError/objectiveNorm_;
2978  if ((1.0-move)*errorCheck>dualTolerance()) {
2979    if (nextGap<gammad*(1.0-move)*errorCheck) {
2980      goodMove=false;
2981    } 
2982  }
2983  if (goodMove)
2984    bestNextGap=nextGap;
2985  return goodMove;
2986}
2987// updateSolution.  Updates solution at end of iteration
2988//returns number fixed
2989int ClpPredictorCorrector::updateSolution(double nextGap)
2990{
2991  int numberTotal = numberRows_+numberColumns_;
2992  //update pi
2993  multiplyAdd(deltaY_,numberRows_,actualDualStep_,dual_,1.0);
2994  CoinZeroN(errorRegion_,numberRows_);
2995  CoinZeroN(rhsFixRegion_,numberRows_);
2996  double maximumRhsInfeasibility=0.0;
2997  double maximumBoundInfeasibility=0.0;
2998  double maximumDualError=1.0e-12;
2999  double primalObjectiveValue=0.0;
3000  double dualObjectiveValue=0.0;
3001  double solutionNorm=1.0e-12;
3002  int numberKilled=0;
3003  double freeMultiplier=1.0e6;
3004  double trueNorm =diagonalNorm_/diagonalScaleFactor_;
3005  if (freeMultiplier<trueNorm) {
3006    freeMultiplier=trueNorm;
3007  } 
3008  if (freeMultiplier>1.0e12) {
3009    freeMultiplier=1.0e12;
3010  } 
3011  freeMultiplier=0.5/freeMultiplier;
3012  double condition = fabs(cholesky_->choleskyCondition());
3013  bool caution;
3014  if ((condition<1.0e10&&trueNorm<1.0e12)||numberIterations_<20) {
3015    caution=false;
3016  } else {
3017    caution=true;
3018  } 
3019  double extra=eExtra;
3020  const double largeFactor=1.0e2;
3021  double largeGap=largeFactor*solutionNorm_;
3022  if (largeGap<largeFactor) {
3023    largeGap=largeFactor;
3024  } 
3025  double dualFake=0.0;
3026  double dualTolerance =  dblParam_[ClpDualTolerance];
3027  dualTolerance=dualTolerance/scaleFactor_;
3028  if (dualTolerance<1.0e-12) {
3029    dualTolerance=1.0e-12;
3030  } 
3031  double offsetObjective=0.0;
3032  const double killTolerance=primalTolerance();
3033  double qDiagonal;
3034  if (mu_<1.0) {
3035    qDiagonal=1.0e-8;
3036  } else {
3037    qDiagonal=1.0e-8*mu_;
3038  } 
3039  //double nextMu = nextGap/(static_cast<double>(2*numberComplementarityPairs_));
3040  //printf("using gap of %g\n",nextMu);
3041  //qDiagonal *= 1.0e2;
3042  //largest allowable ratio of lowerSlack/zVec (etc)
3043  double largestRatio;
3044  double epsilonBase;
3045  double diagonalLimit;
3046  if (!caution) {
3047    epsilonBase=eBase;
3048    largestRatio=eRatio;
3049    diagonalLimit=eDiagonal;
3050  } else {
3051    epsilonBase=eBaseCaution;
3052    largestRatio=eRatioCaution;
3053    diagonalLimit=eDiagonalCaution;
3054  } 
3055  double smallGap=1.0e2;
3056  double maximumDJInfeasibility=0.0;
3057  int numberIncreased=0;
3058  int numberDecreased=0;
3059  double largestDiagonal=0.0;
3060  double smallestDiagonal=1.0e50;
3061  double largeGap2 = CoinMax(1.0e7,1.0e2*solutionNorm_);
3062  //largeGap2 = 1.0e9;
3063  // When to start looking at killing (factor0
3064  double killFactor;
3065#ifndef NO_RTTI
3066  ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
3067#else
3068  ClpQuadraticObjective * quadraticObj = NULL;
3069  if (objective_->type()==2)
3070    quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
3071#endif
3072  if (!quadraticObj||1) {
3073    if (numberIterations_<50) {
3074      killFactor = 1.0;
3075    } else if (numberIterations_<100) {
3076      killFactor = 10.0;
3077      stepLength_=CoinMax(stepLength_,0.9995);
3078    }else if (numberIterations_<150) {
3079      killFactor = 100.0;
3080      stepLength_=CoinMax(stepLength_,0.99995);
3081    } else {
3082      killFactor = 1.0e5;
3083      stepLength_=CoinMax(stepLength_,0.999995);
3084    }
3085  } else {
3086    killFactor=1.0;
3087  }
3088  // put next primal into deltaSL_
3089  int iColumn;
3090  int iRow;
3091  for (iColumn=0;iColumn<numberTotal;iColumn++) {
3092    double thisWeight=deltaX_[iColumn];
3093    double newPrimal=solution_[iColumn]+1.0*actualPrimalStep_*thisWeight;
3094    deltaSL_[iColumn]=newPrimal;
3095  }
3096#if 0
3097  // nice idea but doesn't work
3098  multiplyAdd(solution_+numberColumns_,numberRows_,-1.0,errorRegion_,0.0);
3099  matrix_->times(1.0,solution_,errorRegion_);
3100  multiplyAdd(deltaSL_+numberColumns_,numberRows_,-1.0,rhsFixRegion_,0.0);
3101  matrix_->times(1.0,deltaSL_,rhsFixRegion_);
3102  double newNorm =  maximumAbsElement(deltaSL_,numberTotal);
3103  double tol = newNorm*primalTolerance();
3104  bool goneInf=false;
3105  for (iRow=0;iRow<numberRows_;iRow++) {
3106    double value=errorRegion_[iRow];
3107    double valueNew=rhsFixRegion_[iRow];
3108    if (fabs(value)<tol&&fabs(valueNew)>tol) {
3109      printf("row %d old %g new %g\n",iRow,value,valueNew);
3110      goneInf=true;
3111    }
3112  }
3113  if (goneInf) {
3114    actualPrimalStep_ *= 0.5;
3115    for (iColumn=0;iColumn<numberTotal;iColumn++) {
3116      double thisWeight=deltaX_[iColumn];
3117      double newPrimal=solution_[iColumn]+1.0*actualPrimalStep_*thisWeight;
3118      deltaSL_[iColumn]=newPrimal;
3119    }
3120  }
3121  CoinZeroN(errorRegion_,numberRows_);
3122  CoinZeroN(rhsFixRegion_,numberRows_);
3123#endif
3124  // do reduced costs
3125  CoinMemcpyN(dual_,numberRows_,dj_+numberColumns_);
3126  CoinMemcpyN(cost_,numberColumns_,dj_);
3127  double quadraticOffset=quadraticDjs(dj_,deltaSL_,1.0);
3128  // Save modified costs for fixed variables
3129  CoinMemcpyN(dj_,numberColumns_,deltaSU_);
3130  matrix_->transposeTimes(-1.0,dual_,dj_);
3131  double gamma2 = gamma_*gamma_; // gamma*gamma will be added to diagonal
3132  double gammaOffset=0.0;
3133  for (iColumn=0;iColumn<numberTotal;iColumn++) {
3134    if (!flagged(iColumn)) {
3135      double reducedCost=dj_[iColumn];
3136      bool thisKilled=false;
3137      double zValue = zVec_[iColumn] + actualDualStep_*deltaZ_[iColumn];
3138      double wValue = wVec_[iColumn] + actualDualStep_*deltaW_[iColumn];
3139      zVec_[iColumn]=zValue;
3140      wVec_[iColumn]=wValue;
3141      double thisWeight=deltaX_[iColumn];
3142      double oldPrimal = solution_[iColumn];
3143      double newPrimal=solution_[iColumn]+actualPrimalStep_*thisWeight;
3144      double dualObjectiveThis=0.0;
3145      double sUpper=extra;
3146      double sLower=extra;
3147      double kill;
3148      if (fabs(newPrimal)>1.0e4) {
3149        kill=killTolerance*1.0e-4*newPrimal;
3150      } else {
3151        kill=killTolerance;
3152      } 
3153      kill*=1.0e-3;//be conservative
3154      double smallerSlack=COIN_DBL_MAX;
3155      bool fakeOldBounds=false;
3156      bool fakeNewBounds=false;
3157      double trueLower;
3158      double trueUpper;
3159      if (iColumn<numberColumns_) {
3160        trueLower = columnLower_[iColumn];
3161        trueUpper = columnUpper_[iColumn];
3162      } else {
3163        trueLower = rowLower_[iColumn-numberColumns_];
3164        trueUpper = rowUpper_[iColumn-numberColumns_];
3165      }
3166      if (oldPrimal>trueLower+largeGap2&&
3167          oldPrimal<trueUpper-largeGap2)
3168        fakeOldBounds=true;
3169      if (newPrimal>trueLower+largeGap2&&
3170          newPrimal<trueUpper-largeGap2)
3171        fakeNewBounds=true;
3172      if (fakeOldBounds) {
3173        if (fakeNewBounds) {
3174          lower_[iColumn]=newPrimal-largeGap2;
3175          lowerSlack_[iColumn] = largeGap2;
3176          upper_[iColumn]=newPrimal+largeGap2;
3177          upperSlack_[iColumn] = largeGap2;
3178        } else {
3179          lower_[iColumn]=trueLower;
3180          setLowerBound(iColumn);
3181          lowerSlack_[iColumn] = CoinMax(newPrimal-trueLower,1.0);
3182          upper_[iColumn]=trueUpper;
3183          setUpperBound(iColumn);
3184          upperSlack_[iColumn] = CoinMax(trueUpper-newPrimal,1.0);
3185        }
3186      } else if (fakeNewBounds) {
3187        lower_[iColumn]=newPrimal-largeGap2;
3188        lowerSlack_[iColumn] = largeGap2;
3189        upper_[iColumn]=newPrimal+largeGap2;
3190        upperSlack_[iColumn] = largeGap2;
3191        // so we can just have one test
3192        fakeOldBounds=true;
3193      }
3194      double lowerBoundInfeasibility=0.0;
3195      double upperBoundInfeasibility=0.0;
3196      if (lowerBound(iColumn)) {
3197        double oldSlack = lowerSlack_[iColumn];
3198        double newSlack;
3199        newSlack=
3200          lowerSlack_[iColumn]+actualPrimalStep_*(oldPrimal-oldSlack
3201                                                 + thisWeight-lower_[iColumn]);
3202        if (fakeOldBounds)
3203          newSlack = lowerSlack_[iColumn];
3204        double epsilon = fabs(newSlack)*epsilonBase;
3205        if (epsilon>1.0e-5) {
3206          //cout<<"bad"<<endl;
3207          epsilon=1.0e-5;
3208        } 
3209        //epsilon=1.0e-14;
3210        //make sure reasonable
3211        if (zValue<epsilon) {
3212          zValue=epsilon;
3213        } 
3214        double feasibleSlack=newPrimal-lower_[iColumn];
3215        if (feasibleSlack>0.0&&newSlack>0.0) {
3216          double smallGap2=smallGap;
3217          if (fabs(0.1*newPrimal)>smallGap) {
3218            smallGap2=0.1*fabs(newPrimal);
3219          } 
3220          double larger;
3221          if (newSlack>feasibleSlack) {
3222            larger=newSlack;
3223          } else {
3224            larger=feasibleSlack;
3225          } 
3226          if (fabs(feasibleSlack-newSlack)<1.0e-6*larger) {
3227            newSlack=feasibleSlack;
3228          } 
3229        } 
3230        if (zVec_[iColumn]>dualTolerance) {
3231          dualObjectiveThis+=lower_[iColumn]*zVec_[iColumn];
3232        } 
3233        lowerSlack_[iColumn]=newSlack;
3234        if (newSlack<smallerSlack) {
3235          smallerSlack=newSlack;
3236        } 
3237        lowerBoundInfeasibility = fabs(newPrimal-lowerSlack_[iColumn]-lower_[iColumn]);
3238        if (lowerSlack_[iColumn]<=kill*killFactor&&fabs(newPrimal-lower_[iColumn])<=kill*killFactor) {
3239          double step = CoinMin(actualPrimalStep_*1.1,1.0);
3240          double newPrimal2=solution_[iColumn]+step*thisWeight;
3241          if (newPrimal2<newPrimal&&dj_[iColumn]>1.0e-5&&numberIterations_>50-40) {
3242            newPrimal=lower_[iColumn];
3243            lowerSlack_[iColumn]=0.0;
3244            //printf("fixing %d to lower\n",iColumn);
3245          }
3246        }
3247        if (lowerSlack_[iColumn]<=kill&&fabs(newPrimal-lower_[iColumn])<=kill) {
3248          //may be better to leave at value?
3249          newPrimal=lower_[iColumn];
3250          lowerSlack_[iColumn]=0.0;
3251          thisKilled=true;
3252          //cout<<j<<" l "<<reducedCost<<" "<<zVec_[iColumn]<<endl;
3253        } else {
3254          sLower+=lowerSlack_[iColumn];
3255        } 
3256      } 
3257      if (upperBound(iColumn)) {
3258        double oldSlack = upperSlack_[iColumn];
3259        double newSlack;
3260        newSlack=
3261          upperSlack_[iColumn]+actualPrimalStep_*(-oldPrimal-oldSlack
3262                                                 - thisWeight+upper_[iColumn]);
3263        if (fakeOldBounds)
3264          newSlack = upperSlack_[iColumn];
3265        double epsilon = fabs(newSlack)*epsilonBase;
3266        if (epsilon>1.0e-5) {
3267          //cout<<"bad"<<endl;
3268          epsilon=1.0e-5;
3269        } 
3270        //make sure reasonable
3271        //epsilon=1.0e-14;
3272        if (wValue<epsilon) {
3273          wValue=epsilon;
3274        } 
3275        double feasibleSlack=upper_[iColumn]-newPrimal;
3276        if (feasibleSlack>0.0&&newSlack>0.0) {
3277          double smallGap2=smallGap;
3278          if (fabs(0.1*newPrimal)>smallGap) {
3279            smallGap2=0.1*fabs(newPrimal);
3280          } 
3281          double larger;
3282          if (newSlack>feasibleSlack) {
3283            larger=newSlack;
3284          } else {
3285            larger=feasibleSlack;
3286          } 
3287          if (fabs(feasibleSlack-newSlack)<1.0e-6*larger) {
3288            newSlack=feasibleSlack;
3289          } 
3290        } 
3291        if (wVec_[iColumn]>dualTolerance) {
3292          dualObjectiveThis-=upper_[iColumn]*wVec_[iColumn];
3293        } 
3294        upperSlack_[iColumn]=newSlack;
3295        if (newSlack<smallerSlack) {
3296          smallerSlack=newSlack;
3297        } 
3298        upperBoundInfeasibility = fabs(newPrimal+upperSlack_[iColumn]-upper_[iColumn]);
3299        if (upperSlack_[iColumn]<=kill*killFactor&&fabs(newPrimal-upper_[iColumn])<=kill*killFactor) {
3300          double step = CoinMin(actualPrimalStep_*1.1,1.0);
3301          double newPrimal2=solution_[iColumn]+step*thisWeight;
3302          if (newPrimal2>newPrimal&&dj_[iColumn]<-1.0e-5&&numberIterations_>50-40) {
3303            newPrimal=upper_[iColumn];
3304            upperSlack_[iColumn]=0.0;
3305            //printf("fixing %d to upper\n",iColumn);
3306          }
3307        }
3308        if (upperSlack_[iColumn]<=kill&&fabs(newPrimal-upper_[iColumn])<=kill) {
3309          //may be better to leave at value?
3310          newPrimal=upper_[iColumn];
3311          upperSlack_[iColumn]=0.0;
3312          thisKilled=true;
3313        } else {
3314          sUpper+=upperSlack_[iColumn];
3315        } 
3316      } 
3317      solution_[iColumn]=newPrimal;
3318      if (fabs(newPrimal)>solutionNorm) {
3319        solutionNorm=fabs(newPrimal);
3320      } 
3321      if (!thisKilled) {
3322        double gammaTerm = gamma2;
3323        if (primalR_) {
3324          gammaTerm += primalR_[iColumn];
3325          quadraticOffset += newPrimal*newPrimal*primalR_[iColumn];
3326        }
3327        double dualInfeasibility=
3328          reducedCost-zVec_[iColumn]+wVec_[iColumn]+gammaTerm*newPrimal;
3329        if (fabs(dualInfeasibility)>dualTolerance) {
3330#if 0
3331          if (dualInfeasibility>0.0) {
3332            // To improve we could reduce t and/or increase z
3333            if (lowerBound(iColumn)) {
3334              double complementarity =zVec_[iColumn]*lowerSlack_[iColumn];
3335              if (complementarity<nextMu) {
3336                double change=
3337                  CoinMin(dualInfeasibility,
3338                      (nextMu-complementarity)/lowerSlack_[iColumn]);
3339                dualInfeasibility -= change;
3340                printf("%d lb locomp %g - dual inf from %g to %g\n",
3341                       iColumn,complementarity,dualInfeasibility+change,
3342                       dualInfeasibility);
3343                zVec_[iColumn] += change;
3344                zValue = CoinMax(zVec_[iColumn],1.0e-12);
3345              }
3346            }
3347            if (upperBound(iColumn)) {
3348              double complementarity =wVec_[iColumn]*upperSlack_[iColumn];
3349              if (complementarity>nextMu) {
3350                double change=
3351                  CoinMin(dualInfeasibility,
3352                      (complementarity-nextMu)/upperSlack_[iColumn]);
3353                dualInfeasibility -= change;
3354                printf("%d ub hicomp %g - dual inf from %g to %g\n",
3355                       iColumn,complementarity,dualInfeasibility+change,
3356                       dualInfeasibility);
3357                wVec_[iColumn] -= change;
3358                wValue = CoinMax(wVec_[iColumn],1.0e-12);
3359              }
3360            }
3361          } else {
3362            // To improve we could reduce z and/or increase t
3363            if (lowerBound(iColumn)) {
3364              double complementarity =zVec_[iColumn]*lowerSlack_[iColumn];
3365              if (complementarity>nextMu) {
3366                double change=
3367                  CoinMax(dualInfeasibility,
3368                      (nextMu-complementarity)/lowerSlack_[iColumn]);
3369                dualInfeasibility -= change;
3370                printf("%d lb hicomp %g - dual inf from %g to %g\n",
3371                       iColumn,complementarity,dualInfeasibility+change,
3372                       dualInfeasibility);
3373                zVec_[iColumn] += change;
3374                zValue = CoinMax(zVec_[iColumn],1.0e-12);
3375              }
3376            }
3377            if (upperBound(iColumn)) {
3378              double complementarity =wVec_[iColumn]*upperSlack_[iColumn];
3379              if (complementarity<nextMu) {
3380                double change=
3381                  CoinMax(dualInfeasibility,
3382                      (complementarity-nextMu)/upperSlack_[iColumn]);
3383                dualInfeasibility -= change;
3384                printf("%d ub locomp %g - dual inf from %g to %g\n",
3385                       iColumn,complementarity,dualInfeasibility+change,
3386                       dualInfeasibility);
3387                wVec_[iColumn] -= change;
3388                wValue = CoinMax(wVec_[iColumn],1.0e-12);
3389              }
3390            }
3391          }
3392#endif
3393          dualFake+=newPrimal*dualInfeasibility;
3394        } 
3395        if (lowerBoundInfeasibility>maximumBoundInfeasibility) {
3396          maximumBoundInfeasibility=lowerBoundInfeasibility;
3397        } 
3398        if (upperBoundInfeasibility>maximumBoundInfeasibility) {
3399          maximumBoundInfeasibility=upperBoundInfeasibility;
3400        } 
3401        dualInfeasibility=fabs(dualInfeasibility);
3402        if (dualInfeasibility>maximumDualError) {
3403          //printf("bad dual %d %g\n",iColumn,
3404          // reducedCost-zVec_[iColumn]+wVec_[iColumn]+gammaTerm*newPrimal);
3405          maximumDualError=dualInfeasibility;
3406        } 
3407        dualObjectiveValue+=dualObjectiveThis;
3408        gammaOffset += newPrimal*newPrimal;
3409        if (sLower>largeGap) {
3410          sLower=largeGap;
3411        } 
3412        if (sUpper>largeGap) {
3413          sUpper=largeGap;
3414        } 
3415#if 1
3416        double divisor = sLower*wValue+sUpper*zValue+gammaTerm*sLower*sUpper;
3417        double diagonalValue=(sUpper*sLower)/divisor;
3418#else
3419        double divisor = sLower*wValue+sUpper*zValue+gammaTerm*sLower*sUpper;
3420        double diagonalValue2=(sUpper*sLower)/divisor;
3421        double diagonalValue;
3422        if (!lowerBound(iColumn)) {
3423          diagonalValue = wValue/sUpper + gammaTerm;
3424        } else if (!upperBound(iColumn)) {
3425          diagonalValue = zValue/sLower + gammaTerm;
3426        } else {
3427          diagonalValue = zValue/sLower + wValue/sUpper + gammaTerm;
3428        }
3429        diagonalValue = 1.0/diagonalValue;
3430#endif
3431        diagonal_[iColumn]=diagonalValue;
3432        //FUDGE
3433        if (diagonalValue>diagonalLimit) {
3434#ifdef COIN_DEVELOP
3435          std::cout<<"large diagonal "<<diagonalValue<<std::endl;
3436#endif
3437          diagonal_[iColumn]=diagonalLimit;
3438        } 
3439        if (diagonalValue<1.0e-10) {
3440          //std::cout<<"small diagonal "<<diagonalValue<<std::endl;
3441        } 
3442        if (diagonalValue>largestDiagonal) {
3443          largestDiagonal=diagonalValue;
3444        } 
3445        if (diagonalValue<smallestDiagonal) {
3446          smallestDiagonal=diagonalValue;
3447        } 
3448        deltaX_[iColumn]=0.0;
3449      } else {
3450        numberKilled++;
3451        if (solution_[iColumn]!=lower_[iColumn]&&
3452            solution_[iColumn]!=upper_[iColumn]) {
3453          printf("%d %g %g %g\n",iColumn,lower_[iColumn],
3454                 solution_[iColumn],upper_[iColumn]);
3455        }
3456        diagonal_[iColumn]=0.0;
3457        zVec_[iColumn]=0.0;
3458        wVec_[iColumn]=0.0;
3459        setFlagged(iColumn);
3460        setFixedOrFree(iColumn);
3461        deltaX_[iColumn]=newPrimal;
3462        offsetObjective+=newPrimal*deltaSU_[iColumn];
3463      } 
3464    } else {
3465      deltaX_[iColumn]=solution_[iColumn];
3466      diagonal_[iColumn]=0.0;
3467      offsetObjective+=solution_[iColumn]*deltaSU_[iColumn];
3468      if (upper_[iColumn]-lower_[iColumn]>1.0e-5) {
3469        if (solution_[iColumn]<lower_[iColumn]+1.0e-8&&dj_[iColumn]<-1.0e-8) {
3470          if (-dj_[iColumn]>maximumDJInfeasibility) {
3471            maximumDJInfeasibility=-dj_[iColumn];
3472          } 
3473        } 
3474        if (solution_[iColumn]>upper_[iColumn]-1.0e-8&&dj_[iColumn]>1.0e-8) {
3475          if (dj_[iColumn]>maximumDJInfeasibility) {
3476            maximumDJInfeasibility=dj_[iColumn];
3477          } 
3478        } 
3479      } 
3480    } 
3481    primalObjectiveValue+=solution_[iColumn]*cost_[iColumn];
3482  }
3483  handler_->message(CLP_BARRIER_DIAGONAL,messages_)
3484    <<largestDiagonal<<smallestDiagonal
3485    <<CoinMessageEol;
3486#if 0
3487  // If diagonal wild - kill some
3488  if (largestDiagonal>1.0e17*smallestDiagonal) {
3489    double killValue =largestDiagonal*1.0e-17;
3490    for (int iColumn=0;iColumn<numberTotal;iColumn++) {
3491      if (fabs(diagonal_[iColumn])<killValue)
3492        diagonal_[iolumn]=0.0;
3493    }
3494  }
3495#endif
3496  // update rhs region
3497  multiplyAdd(deltaX_+numberColumns_,numberRows_,-1.0,rhsFixRegion_,1.0);
3498  matrix_->times(1.0,deltaX_,rhsFixRegion_);
3499  primalObjectiveValue += 0.5*gamma2*gammaOffset+0.5*quadraticOffset; 
3500  if (quadraticOffset) {
3501    //  printf("gamma offset %g %g, quadoffset %g\n",gammaOffset,gamma2*gammaOffset,quadraticOffset);
3502  }
3503 
3504  dualObjectiveValue+=offsetObjective+dualFake;
3505  dualObjectiveValue -= 0.5*gamma2*gammaOffset+0.5*quadraticOffset; 
3506  if (numberIncreased||numberDecreased) {
3507    handler_->message(CLP_BARRIER_SLACKS,messages_)
3508      <<numberIncreased<<numberDecreased
3509      <<CoinMessageEol;
3510  } 
3511  if (maximumDJInfeasibility) {
3512    handler_->message(CLP_BARRIER_DUALINF,messages_)
3513      <<maximumDJInfeasibility
3514      <<CoinMessageEol;
3515  } 
3516  // Need to rethink (but it is only for printing)
3517  sumPrimalInfeasibilities_=maximumRhsInfeasibility;
3518  sumDualInfeasibilities_=maximumDualError;
3519  maximumBoundInfeasibility_ = maximumBoundInfeasibility;
3520  //compute error and fixed RHS
3521  multiplyAdd(solution_+numberColumns_,numberRows_,-1.0,errorRegion_,0.0);
3522  matrix_->times(1.0,solution_,errorRegion_);
3523  maximumDualError_=maximumDualError;
3524  maximumBoundInfeasibility_=maximumBoundInfeasibility;
3525  solutionNorm_=solutionNorm;
3526  //finish off objective computation
3527  primalObjective_=primalObjectiveValue*scaleFactor_;
3528  double dualValue2=innerProduct(dual_,numberRows_,
3529                rhsFixRegion_);
3530  dualObjectiveValue-=dualValue2;
3531  dualObjective_=dualObjectiveValue*scaleFactor_;
3532  if (numberKilled) {
3533      handler_->message(CLP_BARRIER_KILLED,messages_)
3534        <<numberKilled
3535        <<CoinMessageEol;
3536  } 
3537  double maximumRHSError1=0.0;
3538  double maximumRHSError2=0.0;
3539  double primalOffset=0.0;
3540  char * dropped = cholesky_->rowsDropped();
3541  for (iRow=0;iRow<numberRows_;iRow++) {
3542    double value=errorRegion_[iRow];
3543    if (!dropped[iRow]) {
3544      if (fabs(value)>maximumRHSError1) {
3545        maximumRHSError1=fabs(value);
3546      } 
3547    } else {
3548      if (fabs(value)>maximumRHSError2) {
3549        maximumRHSError2=fabs(value);
3550      } 
3551      primalOffset+=value*dual_[iRow];
3552    } 
3553  } 
3554  primalObjective_-=primalOffset*scaleFactor_;
3555  if (maximumRHSError1>maximumRHSError2) {
3556    maximumRHSError_=maximumRHSError1;
3557  } else {
3558    maximumRHSError_=maximumRHSError1; //note change
3559    if (maximumRHSError2>primalTolerance()) {
3560      handler_->message(CLP_BARRIER_ABS_DROPPED,messages_)
3561        <<maximumRHSError2
3562        <<CoinMessageEol;
3563    } 
3564    //if (maximumRHSError2>1.0e-1) {
3565      //set all bits false
3566      //cholesky_->resetRowsDropped();
3567    //}
3568  } 
3569  objectiveNorm_=maximumAbsElement(dual_,numberRows_);
3570  if (objectiveNorm_<1.0e-12) {
3571    objectiveNorm_=1.0e-12;
3572  } 
3573  if (objectiveNorm_<baseObjectiveNorm_) {
3574    //std::cout<<" base "<<baseObjectiveNorm_<<" "<<objectiveNorm_<<std::endl;
3575    if (objectiveNorm_<baseObjectiveNorm_*1.0e-4) {
3576      objectiveNorm_=baseObjectiveNorm_*1.0e-4;
3577    } 
3578  } 
3579  bool primalFeasible=true;
3580  if (maximumRHSError_>primalTolerance()||
3581      maximumDualError_>dualTolerance/scaleFactor_) {
3582    handler_->message(CLP_BARRIER_ABS_ERROR,messages_)
3583      <<maximumRHSError_<<maximumDualError_
3584      <<CoinMessageEol;
3585  } 
3586  if (rhsNorm_>solutionNorm_) {
3587    solutionNorm_=rhsNorm_;
3588  } 
3589  double scaledRHSError=maximumRHSError_/(solutionNorm_+10.0);
3590  bool dualFeasible=true;
3591#if KEEP_GOING_IF_FIXED > 5
3592  if (maximumBoundInfeasibility_>primalTolerance()||
3593      scaledRHSError>primalTolerance())
3594    primalFeasible=false;
3595#else
3596  if (maximumBoundInfeasibility_>primalTolerance()||
3597      scaledRHSError>CoinMax(CoinMin(100.0*primalTolerance(),1.0e-5),
3598                             primalTolerance())
3599    primalFeasible=false;
3600#endif
3601  // relax dual test if obj big and gap smallish
3602  double gap=fabs(primalObjective_-dualObjective_);
3603  double sizeObj = CoinMin(fabs(primalObjective_),fabs(dualObjective_))+1.0e-50;
3604  //printf("gap %g sizeObj %g ratio %g comp %g\n",
3605  //     gap,sizeObj,gap/sizeObj,complementarityGap_);
3606  if (numberIterations_>100&&gap/sizeObj<1.0e-9&&complementarityGap_<1.0e-7*sizeObj)
3607    dualTolerance *= 1.0e2;
3608  if (maximumDualError_>objectiveNorm_*dualTolerance) 
3609    dualFeasible=false;
3610  if (!primalFeasible||!dualFeasible) {
3611    handler_->message(CLP_BARRIER_FEASIBLE,messages_)
3612      <<maximumBoundInfeasibility_<<scaledRHSError
3613      <<maximumDualError_/objectiveNorm_
3614      <<CoinMessageEol;
3615  } 
3616  if (!gonePrimalFeasible_) {
3617    gonePrimalFeasible_=primalFeasible;
3618  } else if (!primalFeasible) {
3619    gonePrimalFeasible_=primalFeasible;
3620    if (!numberKilled) {
3621      handler_->message(CLP_BARRIER_GONE_INFEASIBLE,messages_)
3622      <<CoinMessageEol;
3623    } 
3624  } 
3625  if (!goneDualFeasible_) {
3626    goneDualFeasible_=dualFeasible;
3627  } else if (!dualFeasible) {
3628    handler_->message(CLP_BARRIER_GONE_INFEASIBLE,messages_)
3629      <<CoinMessageEol;
3630    goneDualFeasible_=dualFeasible;
3631  } 
3632  //objectiveValue();
3633  if (solutionNorm_>1.0e40) {
3634    std::cout <<"primal off to infinity"<<std::endl;
3635    abort();
3636  } 
3637  if (objectiveNorm_>1.0e40) {
3638    std::cout <<"dual off to infinity"<<std::endl;
3639    abort();
3640  } 
3641  handler_->message(CLP_BARRIER_STEP,messages_)
3642    <<actualPrimalStep_
3643    <<actualDualStep_
3644    <<mu_
3645    <<CoinMessageEol;
3646  numberIterations_++;
3647  return numberKilled;
3648}
3649//  Save info on products of affine deltaSU*deltaW and deltaSL*deltaZ
3650double 
3651ClpPredictorCorrector::affineProduct()
3652{
3653  double product = 0.0;
3654  //IF zVec starts as 0 then deltaZ always zero
3655  //(remember if free then zVec not 0)
3656  //I think free can be done with careful use of boundSlacks to zero
3657  //out all we want
3658  for (int iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
3659    double w3=deltaZ_[iColumn]*deltaX_[iColumn];
3660    double w4=-deltaW_[iColumn]*deltaX_[iColumn];
3661    if (lowerBound(iColumn)) {
3662      w3+=deltaZ_[iColumn]*(solution_[iColumn]-lowerSlack_[iColumn]-lower_[iColumn]);
3663      product+=w3;
3664    } 
3665    if (upperBound(iColumn)) {
3666      w4+=deltaW_[iColumn]*(-solution_[iColumn]-upperSlack_[iColumn]+upper_[iColumn]);
3667      product+=w4;
3668    } 
3669  } 
3670  return product;
3671}
3672//See exactly what would happen given current deltas
3673void 
3674ClpPredictorCorrector::debugMove(int phase,double primalStep, double dualStep)
3675{
3676#ifndef SOME_DEBUG
3677  return;
3678#endif
3679  int numberTotal = numberRows_+numberColumns_;
3680  double * dualNew = ClpCopyOfArray(dual_,numberRows_);
3681  double * errorRegionNew = new double [numberRows_];
3682  double * rhsFixRegionNew = new double [numberRows_];
3683  double * primalNew = ClpCopyOfArray(solution_,numberTotal);
3684  double * djNew = new double[numberTotal];
3685  //update pi
3686  multiplyAdd(deltaY_,numberRows_,dualStep,dualNew,1.0);
3687  // do reduced costs
3688  CoinMemcpyN(dualNew,numberRows_,djNew+numberColumns_);
3689  CoinMemcpyN(cost_,numberColumns_,djNew);
3690  matrix_->transposeTimes(-1.0,dualNew,djNew);
3691  // update x
3692  int iColumn;
3693  for (iColumn=0;iColumn<numberTotal;iColumn++) {
3694    if (!flagged(iColumn)) 
3695      primalNew[iColumn] +=primalStep*deltaX_[iColumn];
3696  }
3697  double quadraticOffset=quadraticDjs(djNew,primalNew,1.0);
3698  CoinZeroN(errorRegionNew,numberRows_);
3699  CoinZeroN(rhsFixRegionNew,numberRows_);
3700  double maximumBoundInfeasibility=0.0;
3701  double maximumDualError=1.0e-12;
3702  double primalObjectiveValue=0.0;
3703  double dualObjectiveValue=0.0;
3704  double solutionNorm=1.0e-12;
3705  const double largeFactor=1.0e2;
3706  double largeGap=largeFactor*solutionNorm_;
3707  if (largeGap<largeFactor) {
3708    largeGap=largeFactor;
3709  } 
3710  double dualFake=0.0;
3711  double dualTolerance =  dblParam_[ClpDualTolerance];
3712  dualTolerance=dualTolerance/scaleFactor_;
3713  if (dualTolerance<1.0e-12) {
3714    dualTolerance=1.0e-12;
3715  }
3716  double newGap=0.0;
3717  double offsetObjective=0.0;
3718  double gamma2 = gamma_*gamma_; // gamma*gamma will be added to diagonal
3719  double gammaOffset=0.0;
3720  double maximumDjInfeasibility=0.0;
3721  for ( iColumn=0;iColumn<numberTotal;iColumn++) {
3722    if (!flagged(iColumn)) {
3723      double reducedCost=djNew[iColumn];
3724      double zValue = zVec_[iColumn] + dualStep*deltaZ_[iColumn];
3725      double wValue = wVec_[iColumn] + dualStep*deltaW_[iColumn];
3726      double thisWeight=deltaX_[iColumn];
3727      double oldPrimal = solution_[iColumn];
3728      double newPrimal=primalNew[iColumn];
3729      double lowerBoundInfeasibility=0.0;
3730      double upperBoundInfeasibility=0.0;
3731      if (lowerBound(iColumn)) {
3732        double oldSlack = lowerSlack_[iColumn];
3733        double newSlack=
3734          lowerSlack_[iColumn]+primalStep*(oldPrimal-oldSlack
3735                                                 + thisWeight-lower_[iColumn]);
3736        if (zValue>dualTolerance) {
3737          dualObjectiveValue+=lower_[iColumn]*zVec_[iColumn];
3738        } 
3739        lowerBoundInfeasibility = fabs(newPrimal-newSlack-lower_[iColumn]);
3740        newGap += newSlack*zValue;
3741      } 
3742      if (upperBound(iColumn)) {
3743        double oldSlack = upperSlack_[iColumn];
3744        double newSlack=
3745          upperSlack_[iColumn]+primalStep*(-oldPrimal-oldSlack
3746                                                 - thisWeight+upper_[iColumn]);
3747        if (wValue>dualTolerance) {
3748          dualObjectiveValue-=upper_[iColumn]*wVec_[iColumn];
3749        } 
3750        upperBoundInfeasibility = fabs(newPrimal+newSlack-upper_[iColumn]);
3751        newGap += newSlack*wValue;
3752      } 
3753      if (fabs(newPrimal)>solutionNorm) {
3754        solutionNorm=fabs(newPrimal);
3755      } 
3756      double gammaTerm = gamma2;
3757      if (primalR_) {
3758        gammaTerm += primalR_[iColumn];
3759        quadraticOffset += newPrimal*newPrimal*primalR_[iColumn];
3760      }
3761      double dualInfeasibility=
3762        reducedCost-zValue+wValue+gammaTerm*newPrimal;
3763      if (fabs(dualInfeasibility)>dualTolerance) {
3764        dualFake+=newPrimal*dualInfeasibility;
3765      } 
3766      if (lowerBoundInfeasibility>maximumBoundInfeasibility) {
3767        maximumBoundInfeasibility=lowerBoundInfeasibility;
3768      } 
3769      if (upperBoundInfeasibility>maximumBoundInfeasibility) {
3770        maximumBoundInfeasibility=upperBoundInfeasibility;
3771      } 
3772      dualInfeasibility=fabs(dualInfeasibility);
3773      if (dualInfeasibility>maximumDualError) {
3774        //printf("bad dual %d %g\n",iColumn,
3775        // reducedCost-zVec_[iColumn]+wVec_[iColumn]+gammaTerm*newPrimal);
3776        maximumDualError=dualInfeasibility;
3777      } 
3778      gammaOffset += newPrimal*newPrimal;
3779      djNew[iColumn]=0.0;
3780    } else {
3781      offsetObjective+=primalNew[iColumn]*cost_[iColumn];
3782      if (upper_[iColumn]-lower_[iColumn]>1.0e-5) {
3783        if (primalNew[iColumn]<lower_[iColumn]+1.0e-8&&djNew[iColumn]<-1.0e-8) {
3784          if (-djNew[iColumn]>maximumDjInfeasibility) {
3785            maximumDjInfeasibility=-djNew[iColumn];
3786          } 
3787        } 
3788        if (primalNew[iColumn]>upper_[iColumn]-1.0e-8&&djNew[iColumn]>1.0e-8) {
3789          if (djNew[iColumn]>maximumDjInfeasibility) {
3790            maximumDjInfeasibility=djNew[iColumn];
3791          } 
3792        } 
3793      } 
3794      djNew[iColumn]=primalNew[iColumn];
3795    } 
3796    primalObjectiveValue+=solution_[iColumn]*cost_[iColumn];
3797  }
3798  // update rhs region
3799  multiplyAdd(djNew+numberColumns_,numberRows_,-1.0,rhsFixRegionNew,1.0);
3800  matrix_->times(1.0,djNew,rhsFixRegionNew);
3801  primalObjectiveValue += 0.5*gamma2*gammaOffset+0.5*quadraticOffset; 
3802  dualObjectiveValue+=offsetObjective+dualFake;
3803  dualObjectiveValue -= 0.5*gamma2*gammaOffset+0.5*quadraticOffset; 
3804  // Need to rethink (but it is only for printing)
3805  //compute error and fixed RHS
3806  multiplyAdd(primalNew+numberColumns_,numberRows_,-1.0,errorRegionNew,0.0);
3807  matrix_->times(1.0,primalNew,errorRegionNew);
3808  //finish off objective computation
3809  double primalObjectiveNew=primalObjectiveValue*scaleFactor_;
3810  double dualValue2=innerProduct(dualNew,numberRows_,
3811                rhsFixRegionNew);
3812  dualObjectiveValue-=dualValue2;
3813  double dualObjectiveNew=dualObjectiveValue*scaleFactor_;
3814  double maximumRHSError1=0.0;
3815  double maximumRHSError2=0.0;
3816  double primalOffset=0.0;
3817  char * dropped = cholesky_->rowsDropped();
3818  int iRow;
3819  for (iRow=0;iRow<numberRows_;iRow++) {
3820    double value=errorRegionNew[iRow];
3821    if (!dropped[iRow]) {
3822      if (fabs(value)>maximumRHSError1) {
3823        maximumRHSError1=fabs(value);
3824      } 
3825    } else {
3826      if (fabs(value)>maximumRHSError2) {
3827        maximumRHSError2=fabs(value);
3828      } 
3829      primalOffset+=value*dualNew[iRow];
3830    } 
3831  } 
3832  primalObjectiveNew-=primalOffset*scaleFactor_;
3833  double maximumRHSError;
3834  if (maximumRHSError1>maximumRHSError2) {
3835    maximumRHSError=maximumRHSError1;
3836  } else {
3837    maximumRHSError=maximumRHSError1; //note change
3838    if (maximumRHSError2>primalTolerance()) {
3839      handler_->message(CLP_BARRIER_ABS_DROPPED,messages_)
3840        <<maximumRHSError2
3841        <<CoinMessageEol;
3842    } 
3843  }
3844  printf("PH %d %g, %g new comp %g, b %g, p %g, d %g\n",phase,
3845         primalStep,dualStep,newGap,maximumBoundInfeasibility,
3846         maximumRHSError,maximumDualError);
3847  if (handler_->logLevel()>1) 
3848    printf("       objs %g %g\n",
3849           primalObjectiveNew,dualObjectiveNew);
3850  if (maximumDjInfeasibility) {
3851    printf(" max dj error on fixed %g\n",
3852           maximumDjInfeasibility);
3853  } 
3854  delete [] dualNew;
3855  delete [] errorRegionNew;
3856  delete [] rhsFixRegionNew;
3857  delete [] primalNew;
3858  delete [] djNew;
3859}
Note: See TracBrowser for help on using the repository browser.