source: trunk/ClpPredictorCorrector.cpp @ 437

Last change on this file since 437 was 425, checked in by forrest, 16 years ago

Improve cholesky

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