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

Last change on this file since 1402 was 1402, checked in by forrest, 10 years ago

get rid of compiler warnings

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