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

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

try and improve stability - also long double interior point

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