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

Last change on this file since 1321 was 1321, checked in by forrest, 13 years ago

out compiler warnings and stability improvements

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