source: branches/BSP/trunk/Clp/src/ClpPredictorCorrector.cpp @ 1069

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

for secondary status on time

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