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

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

changes for cholesky including code from Anshul Gupta

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