source: tags/move-to-subversion/ClpPredictorCorrector.cpp @ 1355

Last change on this file since 1355 was 739, checked in by forrest, 14 years ago

trying to improve

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