source: trunk/ClpPredictorCorrector.cpp @ 369

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

improving interior point code

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 81.5 KB
Line 
1// Copyright (C) 2003, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4/*
5   Implements crude primal dual predictor corrector algorithm
6
7 */
8
9
10#include "CoinPragma.hpp"
11#include <math.h>
12
13#include "CoinHelperFunctions.hpp"
14#include "ClpPredictorCorrector.hpp"
15#include "CoinPackedMatrix.hpp"
16#include "ClpMessage.hpp"
17#include "ClpCholeskyBase.hpp"
18#include "ClpHelperFunctions.hpp"
19#include <cfloat>
20#include <cassert>
21#include <string>
22#include <cstdio>
23#include <iostream>
24
25static double eScale=1.0e57;
26static double eBaseCaution=1.0e-12;
27static double eBase=1.0e-12;
28static double eRatio=1.0e40;
29static double eRatioCaution=1.0e25;
30static double eDiagonal=1.0e25;
31static double eDiagonalCaution=1.0e18;
32static double eExtra=1.0e-12;
33static double eFree =1.0e3;
34
35// main function
36
37int ClpPredictorCorrector::solve ( )
38{
39  problemStatus_=-1;
40  algorithm_=1;
41  //create all regions
42  if (!createWorkingData()) {
43    problemStatus_=4;
44    return 2;
45  }
46  //initializeFeasible(); - this just set fixed flag
47  smallestInfeasibility_=COIN_DBL_MAX;
48  int i;
49  for (i=0;i<LENGTH_HISTORY;i++) 
50    historyInfeasibility_[i]=COIN_DBL_MAX;
51
52  //bool firstTime=true;
53  //firstFactorization(true);
54  if (cholesky_->order(this)) {
55    printf("Not enough memory\n");
56    problemStatus_=4;
57    return -1;
58  }
59  mu_=1.0e10;
60  diagonalScaleFactor_=1.0;
61  //set iterations
62  numberIterations_=-1;
63  int numberTotal = numberRows_+numberColumns_;
64  //initialize solution here
65  if(createSolution()<0) {
66    printf("Not enough memory\n");
67    problemStatus_=4;
68    return -1;
69  }
70  //firstFactorization(false);
71  CoinZeroN(dual_,numberRows_);
72  multiplyAdd(solution_+numberColumns_,numberRows_,-1.0,errorRegion_,0.0);
73  matrix_->times(1.0,solution_,errorRegion_);
74  maximumRHSError_=maximumAbsElement(errorRegion_,numberRows_);
75  maximumBoundInfeasibility_=maximumRHSError_;
76  //double maximumDualError_=COIN_DBL_MAX;
77  //initialize
78  actualDualStep_=0.0;
79  actualPrimalStep_=0.0;
80  gonePrimalFeasible_=false;
81  goneDualFeasible_=false;
82  //bool hadGoodSolution=false;
83  diagonalNorm_=solutionNorm_;
84  mu_=solutionNorm_;
85  int numberFixed=updateSolution();
86  int numberFixedTotal=numberFixed;
87  //int numberRows_DroppedBefore=0;
88  //double extra=eExtra;
89  //double maximumPerturbation=COIN_DBL_MAX;
90  //constants for infeas interior point
91  const double beta2 = 0.99995;
92  const double tau   = 0.00002;
93  double lastComplementarityGap=COIN_DBL_MAX;
94  int lastGoodIteration=0;
95  double bestObjectiveGap=COIN_DBL_MAX;
96  int saveIteration=-1;
97  bool sloppyOptimal=false;
98  double * savePi=NULL;
99  double * savePrimal=NULL;
100  // Extra regions for centering
101  double * saveX = new double[numberTotal];
102  double * saveY = new double[numberRows_];
103  double * saveZ = new double[numberTotal];
104  double * saveT = new double[numberTotal];
105  while (problemStatus_<0) {
106    //#define FULL_DEBUG
107#ifdef FULL_DEBUG
108    {
109      int i;
110      printf("row    pi          artvec       rhsfx\n");
111      for (i=0;i<numberRows_;i++) {
112        printf("%d %g %g %g\n",i,dual_[i],errorRegion_[i],rhsFixRegion_[i]);
113      }
114      printf(" col  dsol  ddiag  dwvec  dzvec dbdslu dbdsll\n");
115      for (i=0;i<numberColumns_+numberRows_;i++) {
116        printf(" %d %g %g %g %g %g %g\n",i,solution_[i],diagonal_[i],tVec_[i],
117               zVec_[i],upperSlack_[i],lowerSlack_[i]);
118      }
119    }
120#endif
121    complementarityGap_=complementarityGap(numberComplementarityPairs_,
122                                           numberComplementarityItems_,0);
123      handler_->message(CLP_BARRIER_ITERATION,messages_)
124    <<numberIterations_
125    <<primalObjective_*optimizationDirection_- dblParam_[ClpObjOffset]
126    << dualObjective_*optimizationDirection_- dblParam_[ClpObjOffset]
127    <<complementarityGap_
128    <<numberFixedTotal
129    <<cholesky_->rank()
130    <<CoinMessageEol;
131    double goodGapChange;
132    if (!sloppyOptimal) {
133      goodGapChange=0.93;
134    } else {
135      goodGapChange=0.7;
136    } 
137    double gapO;
138    double lastGood=bestObjectiveGap;
139    if (gonePrimalFeasible_&&goneDualFeasible_) {
140      double largestObjective;
141      if (fabs(primalObjective_)>fabs(dualObjective_)) {
142        largestObjective = fabs(primalObjective_);
143      } else {
144        largestObjective = fabs(dualObjective_);
145      } 
146      if (largestObjective<1.0) {
147        largestObjective=1.0;
148      } 
149      gapO=fabs(primalObjective_-dualObjective_)/largestObjective;
150      handler_->message(CLP_BARRIER_OBJECTIVE_GAP,messages_)
151        <<gapO
152        <<CoinMessageEol;
153      //start saving best
154      if (gapO<bestObjectiveGap) {
155        saveIteration=numberIterations_;
156        bestObjectiveGap=gapO;
157        if (!savePi) {
158          savePi=new double[numberRows_];
159          savePrimal = new double [numberTotal];
160        } 
161        CoinMemcpyN(dual_,numberRows_,savePi);
162        CoinMemcpyN(solution_,numberTotal,savePrimal);
163      } else if(gapO>2.0*bestObjectiveGap) {
164        //maybe be more sophisticated e.g. re-initialize having
165        //fixed variables and dropped rows
166        //std::cout <<" gap increasing "<<std::endl;
167      } 
168      //std::cout <<"could stop"<<std::endl;
169      //gapO=0.0;
170      if (fabs(primalObjective_-dualObjective_)<dualTolerance()) {
171        gapO=0.0;
172      } 
173    } else {
174      gapO=COIN_DBL_MAX;
175      if (saveIteration>=0) {
176        handler_->message(CLP_BARRIER_GONE_INFEASIBLE,messages_)
177          <<CoinMessageEol;
178        if (sloppyOptimal) {
179          // vaguely optimal
180          double scaledRHSError=maximumRHSError_/solutionNorm_;
181          if (maximumBoundInfeasibility_>1.0e-2||
182              scaledRHSError>1.0e-2||
183              maximumDualError_>objectiveNorm_*1.0e-2) {
184            handler_->message(CLP_BARRIER_EXIT2,messages_)
185              <<saveIteration
186              <<CoinMessageEol;
187            break;
188          }
189        } else {
190          // not close to optimal but check if getting bad
191          double scaledRHSError=maximumRHSError_/solutionNorm_;
192          if (maximumBoundInfeasibility_>1.0e-1||
193              scaledRHSError>1.0e-1||
194              maximumDualError_>objectiveNorm_*1.0e-1) {
195            handler_->message(CLP_BARRIER_EXIT2,messages_)
196              <<saveIteration
197              <<CoinMessageEol;
198            break;
199          }
200        }
201      } 
202    } 
203    if ((gapO<1.0e-6||(gapO<1.0e-4&&complementarityGap_<0.1))&&!sloppyOptimal) {
204      sloppyOptimal=true;
205      handler_->message(CLP_BARRIER_CLOSE_TO_OPTIMAL,messages_)
206        <<numberIterations_<<complementarityGap_
207        <<CoinMessageEol;
208    } 
209    //if (complementarityGap_>=0.98*lastComplementarityGap) {
210    //tryJustPredictor=true;
211    //printf("trying just predictor\n");
212    //}
213    if (complementarityGap_>=1.05*lastComplementarityGap) {
214      handler_->message(CLP_BARRIER_COMPLEMENTARITY,messages_)
215        <<complementarityGap_<<"increasing"
216        <<CoinMessageEol;
217      if (saveIteration>=0&&sloppyOptimal) {
218        handler_->message(CLP_BARRIER_EXIT2,messages_)
219          <<saveIteration
220          <<CoinMessageEol;
221        break;
222      } else {
223        //lastComplementarityGap=complementarityGap_;
224      } 
225    } else if (complementarityGap_<goodGapChange*lastComplementarityGap) {
226      lastGoodIteration=numberIterations_;
227      lastComplementarityGap=complementarityGap_;
228    } else if (numberIterations_-lastGoodIteration>=5&&complementarityGap_<1.0e-3) {
229      handler_->message(CLP_BARRIER_COMPLEMENTARITY,messages_)
230        <<complementarityGap_<<"not decreasing"
231        <<CoinMessageEol;
232      if (gapO>0.75*lastGood) {
233        break;
234      } 
235    } 
236    if (numberIterations_>maximumBarrierIterations_) {
237      handler_->message(CLP_BARRIER_STOPPING,messages_)
238        <<CoinMessageEol;
239      break;
240    } 
241    if (gapO<targetGap_) {
242      problemStatus_=0;
243      handler_->message(CLP_BARRIER_EXIT,messages_)
244        <<" "
245        <<CoinMessageEol;
246        break;//finished
247    } 
248    if (complementarityGap_<1.0e-12) {
249      problemStatus_=0;
250      handler_->message(CLP_BARRIER_EXIT,messages_)
251        <<"- small complementarity gap"
252        <<CoinMessageEol;
253        break;//finished
254    } 
255    if (complementarityGap_<1.0e-10&&gapO<1.0e-10) {
256      problemStatus_=0;
257      handler_->message(CLP_BARRIER_EXIT,messages_)
258        <<"- objective gap and complementarity gap both small"
259        <<CoinMessageEol;
260        break;//finished
261    } 
262    if (gapO<1.0e-9) {
263      double value=gapO*complementarityGap_;
264      value*=actualPrimalStep_;
265      value*=actualDualStep_;
266      //std::cout<<value<<std::endl;
267      if (value<1.0e-17&&numberIterations_>lastGoodIteration) {
268        problemStatus_=0;
269        handler_->message(CLP_BARRIER_EXIT,messages_)
270          <<"- objective gap and complementarity gap both smallish and small steps"
271          <<CoinMessageEol;
272        break;//finished
273      } 
274    } 
275    double nextGap=COIN_DBL_MAX;
276    int nextNumber=0;
277    int nextNumberItems=0;
278    worstDirectionAccuracy_=0.0;
279    int newDropped=0;
280    //Predictor step
281    //prepare for cholesky.  Set up scaled diagonal in deltaX
282    //  ** for efficiency may be better if scale factor known before
283    double norm2=0.0;
284    double maximumValue;
285    getNorms(diagonal_,numberTotal,maximumValue,norm2);
286    diagonalNorm_ = sqrt(norm2/numberComplementarityPairs_);
287    diagonalScaleFactor_=1.0;
288    double maximumAllowable=eScale;
289    //scale so largest is less than allowable ? could do better
290    double factor=0.5;
291    while (maximumValue>maximumAllowable) {
292      diagonalScaleFactor_*=factor;
293      maximumValue*=factor;
294    } /* endwhile */
295    if (diagonalScaleFactor_!=1.0) {
296      handler_->message(CLP_BARRIER_SCALING,messages_)
297        <<"diagonal"<<diagonalScaleFactor_
298        <<CoinMessageEol;
299      diagonalNorm_*=diagonalScaleFactor_;
300    } 
301    multiplyAdd(NULL,numberTotal,0.0,diagonal_,
302                diagonalScaleFactor_);
303    int * rowsDroppedThisTime = new int [numberRows_];
304    newDropped=cholesky_->factorize(diagonal_,rowsDroppedThisTime);
305    if (newDropped>0) {
306      //int newDropped2=cholesky_->factorize(diagonal_,rowsDroppedThisTime);
307      //assert(!newDropped2);
308      if (newDropped<0) {
309        //replace dropped
310        newDropped=-newDropped;
311        //off 1 to allow for reset all
312        newDropped--;
313        //set all bits false
314        cholesky_->resetRowsDropped();
315      } 
316    } else if (newDropped==-1) {
317      printf("Out of memory\n");
318      problemStatus_=4;
319      return -1;
320    } 
321    delete [] rowsDroppedThisTime;
322    if (cholesky_->status()) {
323      std::cout << "bad cholesky?" <<std::endl;
324      abort();
325    }
326    int phase=0; // predictor, corrector , primal dual
327    double directionAccuracy=0.0;
328    bool doCorrector=true;
329    bool goodMove=true;
330    //set up for affine direction
331    setupForSolve(phase);
332    directionAccuracy=findDirectionVector(phase);
333    if (directionAccuracy>worstDirectionAccuracy_) {
334      worstDirectionAccuracy_=directionAccuracy;
335    } 
336    findStepLength(phase);
337    nextGap=complementarityGap(nextNumber,nextNumberItems,1);
338    double affineGap=nextGap;
339    int bestPhase=0;
340    double bestNextGap=nextGap;
341    // ?
342    bestNextGap=max(nextGap,0.8*complementarityGap_);
343    if (complementarityGap_>1.0e-4*numberComplementarityPairs_) {
344      //std::cout <<"predicted duality gap "<<nextGap<<std::endl;
345      double part1=nextGap/numberComplementarityPairs_;
346      double part2=nextGap/complementarityGap_;
347      mu_=part1*part2*part2;
348#if 0
349      double papermu =complementarityGap_/numberComplementarityPairs_;
350      double affmu = nextGap/nextNumber;
351      double sigma = pow(affmu/papermu,3);
352      printf("mu %g, papermu %g, affmu %g, sigma %g sigmamu %g\n",
353             mu_,papermu,affmu,sigma,sigma*papermu);
354#endif 
355      //printf("paper mu %g\n",(nextGap*nextGap*nextGap)/(complementarityGap_*complementarityGap_*
356      //                                            (double) numberComplementarityPairs_));
357    } else {
358      double phi;
359      if (numberComplementarityPairs_<=500) {
360        phi=pow((double) numberComplementarityPairs_,2.0);
361      } else {
362        phi=pow((double) numberComplementarityPairs_,1.5);
363        if (phi<500.0*500.0) {
364          phi=500.0*500.0;
365        } 
366      }
367      mu_=complementarityGap_/phi;
368    } 
369    //save information
370    double product=affineProduct();
371    //#define ALWAYS
372#ifndef ALWAYS
373    //can we do corrector step?
374    double xx= complementarityGap_*(beta2-tau) +product;
375    if (xx>0.0) {
376      double saveMu = mu_;
377      double mu2=numberComplementarityPairs_;
378      mu2=xx/mu2;
379      if (mu2>mu_) {
380        //std::cout<<" could increase to "<<mu2<<std::endl;
381        //was mu2=mu2*0.25;
382        mu2=mu2*0.99;
383        if (mu2<mu_) {
384          mu_=mu2;
385          //std::cout<<" changing mu to "<<mu_<<std::endl;
386        } else {
387          //std::cout<<std::endl;
388        } 
389      } else {
390        //std::cout<<" should decrease to "<<mu2<<std::endl;
391        mu_=0.5*mu2;
392        //std::cout<<" changing mu to "<<mu_<<std::endl;
393      } 
394      handler_->message(CLP_BARRIER_MU,messages_)
395        <<saveMu<<mu_
396        <<CoinMessageEol;
397    } else {
398      //std::cout<<" bad by any standards"<<std::endl;
399    } 
400    if (complementarityGap_*(beta2-tau)+product-mu_*numberComplementarityPairs_<0.0) {
401#ifdef SOME_DEBUG
402      printf("failed 1 product %g mu %g - %g < 0.0, nextGap %g\n",product,mu_,
403             complementarityGap_*(beta2-tau)+product-mu_*numberComplementarityPairs_,
404             nextGap);
405#endif
406      doCorrector=false;
407      if (nextGap>0.9*complementarityGap_||1) {
408        goodMove=false;
409        bestNextGap=COIN_DBL_MAX;
410      }
411      //double floatNumber = 2.0*numberComplementarityPairs_;
412      //floatNumber = 1.0*numberComplementarityItems_;
413      //mu_=nextGap/floatNumber;
414      handler_->message(CLP_BARRIER_INFO,messages_)
415        <<"no corrector step"
416        <<CoinMessageEol;
417    } else {
418      phase=1;
419    }
420#else
421    phase=1;
422#endif
423    if (goodMove&&doCorrector) {
424      //set up for next step
425      setupForSolve(phase);
426      double directionAccuracy2=findDirectionVector(phase);
427      if (directionAccuracy2>worstDirectionAccuracy_) {
428        worstDirectionAccuracy_=directionAccuracy2;
429      } 
430      double testValue=1.0e2*directionAccuracy;
431      if (1.0e2*projectionTolerance_>testValue) {
432        testValue=1.0e2*projectionTolerance_;
433      } 
434      if (primalTolerance()>testValue) {
435        testValue=primalTolerance();
436      } 
437      if (maximumRHSError_>testValue) {
438        testValue=maximumRHSError_;
439      } 
440      if (directionAccuracy2>testValue&&numberIterations_>=-77) {
441        goodMove=false;
442#ifdef SOME_DEBUG
443        printf("accuracy %g phase 1 failed, test value %g\n",
444               directionAccuracy2,testValue);
445#endif
446        bestNextGap=COIN_DBL_MAX;
447      } 
448      if (goodMove) {
449        findStepLength(phase);
450        // just for debug
451        nextGap = complementarityGap(nextNumber,nextNumberItems,1);
452#ifndef ALWAYS
453        if (numberIterations_>=-77) {
454          goodMove=checkGoodMove(true,bestNextGap);
455          if (!goodMove) {
456#ifdef SOME_DEBUG
457            printf("checkGoodMove failed\n");
458#endif
459            if (affineGap<0.5*complementarityGap_) {
460              //||(affineGap<0.95*complementarityGap_&&(numberIterations_%3)==0)) {
461              // Back to affine
462              setupForSolve(0);
463              directionAccuracy=findDirectionVector(0);
464              findStepLength(0);
465              nextGap=complementarityGap(nextNumber,nextNumberItems,1);
466              goodMove=true;
467            }
468          }
469        } else {
470          goodMove=true;
471        } 
472#endif
473      }
474    } 
475    if (!goodMove) {
476      // Just primal dual step
477      double floatNumber;
478      floatNumber = 2.0*numberComplementarityPairs_;
479      //floatNumber = 1.1*numberComplementarityPairs_;
480      mu_=complementarityGap_/floatNumber;
481      //set up for next step
482      setupForSolve(2);
483      double directionAccuracy=findDirectionVector(2);
484      if (directionAccuracy>worstDirectionAccuracy_) {
485        worstDirectionAccuracy_=directionAccuracy;
486      } 
487      double testValue=1.0e2*directionAccuracy;
488      if (1.0e2*projectionTolerance_>testValue) {
489        testValue=1.0e2*projectionTolerance_;
490      } 
491      if (primalTolerance()>testValue) {
492        testValue=primalTolerance();
493      } 
494      if (maximumRHSError_>testValue) {
495        testValue=maximumRHSError_;
496      } 
497      findStepLength(2);
498      // just for debug
499      nextGap=complementarityGap(nextNumber,nextNumberItems,2);
500      if (nextGap>bestNextGap&&bestPhase==0&&0) {
501        // Back to affine
502        setupForSolve(0);
503        directionAccuracy=findDirectionVector(0);
504        findStepLength(0);
505        nextGap=complementarityGap(nextNumber,nextNumberItems,1);
506      }
507    }
508    if (!goodMove)
509      mu_=nextGap / ((double) 1.1*nextNumber);
510    goodMove=true;
511    // Do centering steps
512    int numberTries=0;
513    double nextCenterGap=0.0;
514    int numberGoodTries=0;
515    double originalDualStep=actualDualStep_;
516    double originalPrimalStep=actualPrimalStep_;
517    while (goodMove&&numberTries<4) {
518      goodMove=false;
519      numberTries++;
520      memcpy(saveX,deltaX_,numberTotal*sizeof(double));
521      memcpy(saveY,deltaY_,numberRows_*sizeof(double));
522      memcpy(saveZ,deltaZ_,numberTotal*sizeof(double));
523      memcpy(saveT,deltaT_,numberTotal*sizeof(double));
524      double savePrimalStep = actualPrimalStep_;
525      double saveDualStep = actualDualStep_;
526      double saveMu = mu_;
527      setupForSolve(3);
528      findDirectionVector(3);
529      findStepLength(3);
530      double xGap = complementarityGap(nextNumber,nextNumberItems,3);
531      // If one small then that's the one that counts
532      double checkDual=saveDualStep;
533      double checkPrimal=savePrimalStep;
534      if (checkDual>5.0*checkPrimal) {
535        checkDual=2.0*checkPrimal;
536      } else if (checkPrimal>5.0*checkDual) {
537        checkPrimal=2.0*checkDual;
538      }
539      //if (actualPrimalStep_<checkPrimal||
540      //  actualDualStep_<checkDual||
541      //  (xGap>nextGap&&xGap>0.9*complementarityGap_)) {
542      if (actualPrimalStep_<checkPrimal||
543          actualDualStep_<checkDual) {
544#ifdef SOME_DEBUG
545        printf("PP rejected gap %g, steps %g %g, 2 gap %g, steps %g %g\n",xGap,
546               actualPrimalStep_,actualDualStep_,nextGap,savePrimalStep,saveDualStep);
547#endif
548        mu_=saveMu;
549        actualPrimalStep_ = savePrimalStep;
550        actualDualStep_ = saveDualStep;
551        memcpy(deltaX_,saveX,numberTotal*sizeof(double));
552        memcpy(deltaY_,saveY,numberRows_*sizeof(double));
553        memcpy(deltaZ_,saveZ,numberTotal*sizeof(double));
554        memcpy(deltaT_,saveT,numberTotal*sizeof(double));
555      } else {
556#ifdef SOME_DEBUG
557        printf("PPphase 3 gap %g, steps %g %g, 2 gap %g, steps %g %g\n",xGap,
558               actualPrimalStep_,actualDualStep_,nextGap,savePrimalStep,saveDualStep);
559#endif
560        numberGoodTries++;
561        nextCenterGap=xGap;
562        // See if big enough change
563        if (actualPrimalStep_<1.01*checkPrimal||
564            actualDualStep_<1.01*checkDual) {
565          // stop now
566        } else {
567          // carry on
568          goodMove=true;
569        }
570      }
571    }
572    if (numberGoodTries&&handler_->logLevel()>1) {
573      printf("%d centering steps moved from (gap %g, dual %g, primal %g) to (gap %g, dual %g, primal %g)\n",
574             numberGoodTries,nextGap,originalDualStep,originalPrimalStep,
575             nextCenterGap, actualDualStep_,actualPrimalStep_);
576    }
577    numberFixed=updateSolution();
578    numberFixedTotal+=numberFixed;
579  } /* endwhile */
580  delete [] saveX;
581  delete [] saveY;
582  delete [] saveZ;
583  delete [] saveT;
584  if (savePi) {
585    //std::cout<<"Restoring from iteration "<<saveIteration<<std::endl;
586    CoinMemcpyN(savePi,numberRows_,dual_);
587    CoinMemcpyN(savePrimal,numberTotal,solution_);
588    delete [] savePi;
589    delete [] savePrimal;
590  } 
591  //recompute slacks
592  // Split out solution
593  CoinZeroN(rowActivity_,numberRows_);
594  CoinMemcpyN(solution_,numberColumns_,columnActivity_);
595  matrix_->times(1.0,columnActivity_,rowActivity_);
596  //unscale objective
597  multiplyAdd(NULL,numberTotal,0.0,cost_,scaleFactor_);
598  multiplyAdd(NULL,numberRows_,0,dual_,scaleFactor_);
599  CoinMemcpyN(cost_,numberColumns_,reducedCost_);
600  matrix_->transposeTimes(-1.0,dual_,reducedCost_);
601  CoinMemcpyN(reducedCost_,numberColumns_,dj_);
602  checkSolution();
603  handler_->message(CLP_BARRIER_END,messages_)
604    <<sumPrimalInfeasibilities_
605    <<sumDualInfeasibilities_
606    <<complementarityGap_
607    <<objectiveValue()
608    <<CoinMessageEol;
609  //#ifdef SOME_DEBUG
610  if (handler_->logLevel()>1) 
611    printf("ENDRUN status %d after %d iterations\n",problemStatus_,numberIterations_);
612  //#endif
613  //std::cout<<"Absolute primal infeasibility at end "<<sumPrimalInfeasibilities_<<std::endl;
614  //std::cout<<"Absolute dual infeasibility at end "<<sumDualInfeasibilities_<<std::endl;
615  //std::cout<<"Absolute complementarity at end "<<complementarityGap_<<std::endl;
616  //std::cout<<"Primal objective "<<objectiveValue()<<std::endl;
617  //std::cout<<"maximum complementarity "<<worstComplementarity_<<std::endl;
618  //delete all temporary regions
619  deleteWorkingData();
620  return problemStatus_;
621}
622// findStepLength.
623//phase  - 0 predictor
624//         1 corrector
625//         2 primal dual
626double ClpPredictorCorrector::findStepLength( int phase)
627{
628  double directionNorm=0.0;
629  double maximumPrimalStep=COIN_DBL_MAX;
630  double maximumDualStep=COIN_DBL_MAX;
631  int numberTotal = numberRows_+numberColumns_;
632  double tolerance = 1.0e-12;
633  int chosenPrimalSequence=-1;
634  int chosenDualSequence=-1;
635  double * zVec = zVec_;
636  double * tVec = tVec_;
637  double * primal = solution_;
638  double * lowerSlack = lowerSlack_;
639  double * upperSlack = upperSlack_;
640  //direction vector in deltaX
641  double * deltaX = deltaX_;
642  int iColumn;
643  for (iColumn=0;iColumn<numberTotal;iColumn++) {
644    if (!flagged(iColumn)) {
645      double z1=0.0;
646      double t1=0.0;
647      double directionElement=deltaX[iColumn];
648      double value=primal[iColumn];
649      if (directionNorm<fabs(directionElement)) {
650        directionNorm=fabs(directionElement);
651      } 
652      if (lowerBound(iColumn)||
653          upperBound(iColumn)) {
654        if (lowerBound(iColumn)) {
655          double delta = - deltaSL_[iColumn];
656          z1 = deltaZ_[iColumn];
657          if (lowerSlack[iColumn]<maximumPrimalStep*delta) {
658            maximumPrimalStep=lowerSlack[iColumn]/delta;
659            chosenPrimalSequence=iColumn;
660          } 
661          if (zVec[iColumn]>tolerance) {
662            if (zVec[iColumn]<-z1*maximumDualStep) {
663              maximumDualStep=-zVec[iColumn]/z1;
664              chosenDualSequence=iColumn;
665            } 
666          } 
667        } else {
668          assert (!zVec[iColumn]);
669        }
670        if (upperBound(iColumn)) {
671          double delta = - deltaSU_[iColumn];;
672          t1 = deltaT_[iColumn];
673          if (upperSlack[iColumn]<maximumPrimalStep*delta) {
674            maximumPrimalStep=upperSlack[iColumn]/delta;
675            chosenPrimalSequence=iColumn;
676          } 
677          if (tVec[iColumn]>tolerance) {
678            if (tVec[iColumn]<-t1*maximumDualStep) {
679              maximumDualStep=-tVec[iColumn]/t1;
680              chosenDualSequence=iColumn;
681            } 
682          } 
683        } else {
684          assert (!tVec[iColumn]);
685        } 
686      } else {
687        //free
688        double gap=fabs(value);
689        double multiplier=1.0/gap;
690        if (gap<1.0) {
691          multiplier=1,0;
692        } 
693        z1 = -zVec[iColumn]-multiplier*directionElement*zVec[iColumn];
694        t1 = -tVec[iColumn] + multiplier*directionElement*tVec[iColumn];
695        deltaZ_[iColumn]=z1;
696        deltaT_[iColumn]=t1;
697      } 
698    } 
699  }
700#ifdef SOME_DEBUG
701  printf("new step - phase %d, norm %g, dual step %g, primal step %g\n",
702         phase,directionNorm,maximumPrimalStep,maximumDualStep);
703#endif
704  actualPrimalStep_=stepLength_*maximumPrimalStep;
705  if (phase>=0&&actualPrimalStep_>1.0) {
706    actualPrimalStep_=1.0;
707  } 
708  actualDualStep_=stepLength_*maximumDualStep;
709  if (phase>=0&&actualDualStep_>1.0) {
710    actualDualStep_=1.0;
711  } 
712#ifdef FULL_DEBUG
713  if (phase==3){
714    double minBeta = 0.1*mu_;
715    double maxBeta = 10.0*mu_;
716    for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
717      if (!flagged(iColumn)) {
718        if (lowerBound(iColumn)) {
719          double change = rhsL_[iColumn] + deltaX_[iColumn];
720          double dualValue=zVec[iColumn]+actualDualStep_*deltaZ_[iColumn];
721          double primalValue=lowerSlack[iColumn]+actualPrimalStep_*change;
722          double gapProduct=dualValue*primalValue;
723          if (delta2Z_[iColumn]<minBeta||delta2Z_[iColumn]>maxBeta)
724            printf("3lower %d primal %g, dual %g, gap %g, old gap %g\n",
725                   iColumn,primalValue,dualValue,gapProduct,delta2Z_[iColumn]);
726        } 
727        if (upperBound(iColumn)) {
728          double change = rhsU_[iColumn]-deltaX_[iColumn];
729          double dualValue=tVec[iColumn]+actualDualStep_*deltaT_[iColumn];
730          double primalValue=upperSlack[iColumn]+actualPrimalStep_*change;
731          double gapProduct=dualValue*primalValue;
732          if (delta2T_[iColumn]<minBeta||delta2T_[iColumn]>maxBeta)
733            printf("3upper %d primal %g, dual %g, gap %g, old gap %g\n",
734                 iColumn,primalValue,dualValue,gapProduct,delta2T_[iColumn]);
735        } 
736      } 
737    }
738  }
739#endif
740#ifdef SOME_DEBUG
741  {
742    double largestL=0.0;
743    double smallestL=COIN_DBL_MAX;
744    double largestU=0.0;
745    double smallestU=COIN_DBL_MAX;
746    double sumL=0.0;
747    double sumU=0.0;
748    int nL=0;
749    int nU=0;
750    for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
751      if (!flagged(iColumn)) {
752        if (lowerBound(iColumn)) {
753          double change = rhsL_[iColumn] + deltaX_[iColumn];
754          double dualValue=zVec[iColumn]+actualDualStep_*deltaZ_[iColumn];
755          double primalValue=lowerSlack[iColumn]+actualPrimalStep_*change;
756          double gapProduct=dualValue*primalValue;
757          largestL = max(largestL,gapProduct);
758          smallestL = min(smallestL,gapProduct);
759          nL++;
760          sumL += gapProduct;
761        } 
762        if (upperBound(iColumn)) {
763          double change = rhsU_[iColumn]-deltaX_[iColumn];
764          double dualValue=tVec[iColumn]+actualDualStep_*deltaT_[iColumn];
765          double primalValue=upperSlack[iColumn]+actualPrimalStep_*change;
766          double gapProduct=dualValue*primalValue;
767          largestU = max(largestU,gapProduct);
768          smallestU = min(smallestU,gapProduct);
769          nU++;
770          sumU += gapProduct;
771        } 
772      } 
773    }
774    double mu = (sumL+sumU)/((double) (nL+nU));
775
776    double minBeta = 0.1*mu;
777    double maxBeta = 10.0*mu;
778    int nBL=0;
779    int nAL=0;
780    int nBU=0;
781    int nAU=0;
782    for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
783      if (!flagged(iColumn)) {
784        if (lowerBound(iColumn)) {
785          double change = rhsL_[iColumn] + deltaX_[iColumn];
786          double dualValue=zVec[iColumn]+actualDualStep_*deltaZ_[iColumn];
787          double primalValue=lowerSlack[iColumn]+actualPrimalStep_*change;
788          double gapProduct=dualValue*primalValue;
789          if (gapProduct<minBeta)
790            nBL++;
791          else if (gapProduct>maxBeta)
792            nAL++;
793        } 
794        if (upperBound(iColumn)) {
795          double change = rhsU_[iColumn]-deltaX_[iColumn];
796          double dualValue=tVec[iColumn]+actualDualStep_*deltaT_[iColumn];
797          double primalValue=upperSlack[iColumn]+actualPrimalStep_*change;
798          double gapProduct=dualValue*primalValue;
799          if (gapProduct<minBeta)
800            nBU++;
801          else if (gapProduct>maxBeta)
802            nAU++;
803        } 
804      } 
805    }
806    printf("phase %d new mu %g new gap %g\n",phase,mu,sumL+sumU);
807    printf("          %d lower, smallest %g, %d below - largest %g, %d above\n",
808           nL,smallestL,nBL,largestL,nAL);
809    printf("          %d upper, smallest %g, %d below - largest %g, %d above\n",
810           nU,smallestU,nBU,largestU,nAU);
811  }
812#endif
813  return directionNorm;
814}
815//#define KKT 2
816/* Does solve. region1 is for deltaX (columns+rows), region2 for deltaPi (rows) */
817void 
818ClpPredictorCorrector::solveSystem(double * region1, double * region2,
819                                   const double * region1In, const double * region2In,
820                                   const double * saveRegion1, const double * saveRegion2,
821                                   bool gentleRefine)
822{
823  int iRow;
824  int numberTotal = numberRows_+numberColumns_;
825  if (region2In) {
826    // normal
827    for (iRow=0;iRow<numberRows_;iRow++)
828      region2[iRow] = region2In[iRow];
829  } else {
830    // initial solution - (diagonal is 1 or 0)
831    CoinZeroN(region2,numberRows_);
832  }
833  int iColumn;
834#if KKT<2
835  for (iColumn=0;iColumn<numberTotal;iColumn++)
836    region1[iColumn] = region1In[iColumn]*diagonal_[iColumn];
837  multiplyAdd(region1+numberColumns_,numberRows_,-1.0,region2,1.0);
838  matrix_->times(1.0,region1,region2);
839  double maximumRHS = maximumAbsElement(region2,numberRows_);
840  double scale=1.0;
841  double unscale=1.0;
842  if (maximumRHS>1.0e-30) {
843    if (maximumRHS<=0.5) {
844      double factor=2.0;
845      while (maximumRHS<=0.5) {
846        maximumRHS*=factor;
847        scale*=factor;
848      } /* endwhile */
849    } else if (maximumRHS>=2.0&&maximumRHS<=COIN_DBL_MAX) {
850      double factor=0.5;
851      while (maximumRHS>=2.0) {
852        maximumRHS*=factor;
853        scale*=factor;
854      } /* endwhile */
855    } 
856    unscale=diagonalScaleFactor_/scale;
857  } else {
858    //effectively zero
859    scale=0.0;
860    unscale=0.0;
861  } 
862  multiplyAdd(NULL,numberRows_,0.0,region2,scale);
863  cholesky_->solve(region2);
864  multiplyAdd(NULL,numberRows_,0.0,region2,unscale);
865  multiplyAdd(region2,numberRows_,-1.0,region1+numberColumns_,0.0);
866  CoinZeroN(region1,numberColumns_);
867  matrix_->transposeTimes(1.0,region2,region1);
868  for (iColumn=0;iColumn<numberTotal;iColumn++)
869    region1[iColumn] = (region1[iColumn]-region1In[iColumn])*diagonal_[iColumn];
870#else
871  for (iColumn=0;iColumn<numberTotal;iColumn++)
872    region1[iColumn] = region1In[iColumn];
873  cholesky_->solveKKT(region1,region2,diagonal_,diagonalScaleFactor_);
874#endif
875  if (saveRegion2) {
876    //refine?
877    double scaleX=1.0;
878    if (gentleRefine) 
879      scaleX=0.8;
880    multiplyAdd(saveRegion2,numberRows_,1.0,region2,scaleX);
881    assert (saveRegion1);
882    multiplyAdd(saveRegion1,numberTotal,1.0,region1,scaleX);
883  } 
884}
885// findDirectionVector.
886double ClpPredictorCorrector::findDirectionVector(const int phase)
887{
888  double projectionTolerance=projectionTolerance_;
889  //temporary
890  //projectionTolerance=1.0e-15;
891  double errorCheck=0.9*maximumRHSError_/solutionNorm_;
892  if (errorCheck>primalTolerance()) {
893    if (errorCheck<projectionTolerance) {
894      projectionTolerance=errorCheck;
895    } 
896  } else {
897    if (primalTolerance()<projectionTolerance) {
898      projectionTolerance=primalTolerance();
899    } 
900  } 
901  double * newError = new double [numberRows_];
902  double * workArray = workArray_;
903  int numberTotal = numberRows_+numberColumns_;
904  //if flagged then entries zero so can do
905  // For KKT separate out
906#ifndef KKT
907  int iColumn;
908  for (iColumn=0;iColumn<numberTotal;iColumn++)
909    deltaX_[iColumn] = workArray[iColumn] - solution_[iColumn];
910  multiplyAdd(deltaX_+numberColumns_,numberRows_,-1.0,deltaY_,0.0);
911  matrix_->times(1.0,deltaX_,deltaY_);
912#else
913  // regions in will be workArray and newError
914  // regions out deltaX_ and deltaY_
915  multiplyAdd(solution_+numberColumns_,numberRows_,1.0,newError,0.0);
916  matrix_->times(-1.0,solution_,newError);
917  double * region1Save=NULL;//for refinement
918#endif
919  bool goodSolve=false;
920  double * regionSave=NULL;//for refinement
921  int numberTries=0;
922  double relativeError=COIN_DBL_MAX;
923  double tryError=1.0e31;
924  while (!goodSolve&&numberTries<30) {
925    double lastError=relativeError;
926    goodSolve=true;
927#ifndef KKT
928    double maximumRHS = maximumAbsElement(deltaY_,numberRows_);
929    double saveMaximum = maximumRHS;
930    double scale=1.0;
931    double unscale=1.0;
932    if (maximumRHS>1.0e-30) {
933      if (maximumRHS<=0.5) {
934        double factor=2.0;
935        while (maximumRHS<=0.5) {
936          maximumRHS*=factor;
937          scale*=factor;
938        } /* endwhile */
939      } else if (maximumRHS>=2.0&&maximumRHS<=COIN_DBL_MAX) {
940        double factor=0.5;
941        while (maximumRHS>=2.0) {
942          maximumRHS*=factor;
943          scale*=factor;
944        } /* endwhile */
945      } 
946      unscale=diagonalScaleFactor_/scale;
947    } else {
948      //effectively zero
949      scale=0.0;
950      unscale=0.0;
951    } 
952    //printf("--putting scales to 1.0\n");
953    //scale=1.0;
954    //unscale=1.0;
955    multiplyAdd(NULL,numberRows_,0.0,deltaY_,scale);
956    cholesky_->solve(deltaY_);
957    multiplyAdd(NULL,numberRows_,0.0,deltaY_,unscale);
958    if (numberTries) {
959      //refine?
960      double scaleX=1.0;
961      if (lastError>1.0e-5) 
962        scaleX=0.8;
963      multiplyAdd(regionSave,numberRows_,1.0,deltaY_,scaleX);
964    } 
965    //CoinZeroN(newError,numberRows_);
966    multiplyAdd(deltaY_,numberRows_,-1.0,deltaX_+numberColumns_,0.0);
967    CoinZeroN(deltaX_,numberColumns_);
968    matrix_->transposeTimes(1.0,deltaY_,deltaX_);
969    //if flagged then entries zero so can do
970    for (iColumn=0;iColumn<numberTotal;iColumn++)
971      deltaX_[iColumn] = deltaX_[iColumn]*diagonal_[iColumn]
972        -workArray[iColumn];
973#else
974    solveSystem(deltaX_, deltaY_,
975                workArray,newError,region1Save,regionSave,lastError>1.0e-5);
976#endif
977    multiplyAdd(deltaX_+numberColumns_,numberRows_,-1.0,newError,0.0);
978    matrix_->times(1.0,deltaX_,newError);
979    numberTries++;
980   
981    //now add in old Ax - doing extra checking
982    double maximumRHSError=0.0;
983    double maximumRHSChange=0.0;
984    int iRow;
985    char * dropped = cholesky_->rowsDropped();
986    for (iRow=0;iRow<numberRows_;iRow++) {
987      if (!dropped[iRow]) {
988        double newValue=newError[iRow];
989        double oldValue=errorRegion_[iRow];
990        //severity of errors depend on signs
991        //**later                                                             */
992        if (fabs(newValue)>maximumRHSChange) {
993          maximumRHSChange=fabs(newValue);
994        } 
995        double result=newValue+oldValue;
996        if (fabs(result)>maximumRHSError) {
997          maximumRHSError=fabs(result);
998        } 
999        newError[iRow]=result;
1000      } else {
1001        double newValue=newError[iRow];
1002        double oldValue=errorRegion_[iRow];
1003        if (fabs(newValue)>maximumRHSChange) {
1004          maximumRHSChange=fabs(newValue);
1005        } 
1006        double result=newValue+oldValue;
1007        newError[iRow]=result;
1008        //newError[iRow]=0.0;
1009        //assert(deltaY_[iRow]==0.0);
1010        deltaY_[iRow]=0.0;
1011      } 
1012    } 
1013    relativeError = maximumRHSError/solutionNorm_;
1014    relativeError = maximumRHSError/saveMaximum;
1015    if (relativeError>tryError) 
1016      relativeError=tryError;
1017    if (relativeError<lastError) {
1018      maximumRHSChange_= maximumRHSChange;
1019      if (relativeError>1.0e-9
1020          ||numberTries>1) {
1021        handler_->message(CLP_BARRIER_ACCURACY,messages_)
1022          <<phase<<numberTries<<relativeError
1023          <<CoinMessageEol;
1024      } 
1025      if (relativeError>projectionTolerance&&numberTries<=3) {
1026        //try and refine
1027        goodSolve=false;
1028      } 
1029      //*** extra test here
1030      if (!goodSolve) {
1031        if (!regionSave) {
1032          regionSave = new double [numberRows_];
1033#ifdef KKT
1034          region1Save = new double [numberTotal];
1035#endif
1036        } 
1037        CoinMemcpyN(deltaY_,numberRows_,regionSave);
1038#ifndef KKT
1039        multiplyAdd(newError,numberRows_,-1.0,deltaY_,0.0);
1040#else
1041        CoinMemcpyN(deltaX_,numberTotal,region1Save);
1042        // and back to input region
1043        CoinMemcpyN(deltaY_,numberRows_,newError);
1044#endif
1045      } 
1046    } else {
1047      //std::cout <<" worse residual = "<<relativeError;
1048      //bring back previous
1049      relativeError=lastError;
1050      CoinMemcpyN(regionSave,numberRows_,deltaY_);
1051#ifndef KKT
1052      multiplyAdd(deltaY_,numberRows_,-1.0,deltaX_+numberColumns_,0.0);
1053      CoinZeroN(deltaX_,numberColumns_);
1054      matrix_->transposeTimes(1.0,deltaY_,deltaX_);
1055      //if flagged then entries zero so can do
1056      for (iColumn=0;iColumn<numberTotal;iColumn++)
1057        deltaX_[iColumn] = deltaX_[iColumn]*diagonal_[iColumn]
1058          -workArray[iColumn];
1059#else
1060        CoinMemcpyN(region1Save,numberTotal,deltaX_);
1061#endif
1062    } 
1063  } /* endwhile */
1064  delete [] regionSave;
1065#ifdef KKT
1066  delete [] region1Save;
1067#endif
1068  delete [] newError;
1069  // now rest
1070  double * zVec = zVec_;
1071  double * tVec = tVec_;
1072  double * lowerSlack = lowerSlack_;
1073  double * upperSlack = upperSlack_;
1074  double extra=eExtra;
1075
1076  for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
1077    deltaSU_[iColumn]=0.0;
1078    deltaSL_[iColumn]=0.0;
1079    deltaZ_[iColumn]=0.0;
1080    deltaT_[iColumn]=0.0;
1081    if (!flagged(iColumn)) {
1082      double deltaX = deltaX_[iColumn];
1083      if (lowerBound(iColumn)||upperBound(iColumn)) {
1084        if (lowerBound(iColumn)) {
1085          double deltaSL = rhsL_[iColumn]+deltaX;
1086          double slack = lowerSlack[iColumn]+extra;
1087          deltaSL_[iColumn]=deltaSL;
1088          deltaZ_[iColumn]=(rhsZ_[iColumn]-zVec[iColumn]*deltaSL)/slack;
1089        } 
1090        if (upperBound(iColumn)) {
1091          double deltaSU = rhsU_[iColumn]-deltaX;
1092          double slack = upperSlack[iColumn]+extra;
1093          deltaSU_[iColumn]=deltaSU;
1094          deltaT_[iColumn]=(rhsT_[iColumn]-tVec[iColumn]*deltaSU)/slack;
1095        }
1096      } else {
1097        // free
1098        double gap=fabs(solution_[iColumn]);
1099        double multiplier=1.0/gap;
1100        if (gap<1.0) {
1101          multiplier=1,0;
1102        } 
1103        deltaZ_[iColumn] = -zVec[iColumn] - multiplier*deltaX*zVec[iColumn];
1104        deltaT_[iColumn] = -tVec[iColumn] + multiplier*deltaX*tVec[iColumn];
1105      }
1106    }
1107  } 
1108  return relativeError;
1109}
1110// createSolution.  Creates solution from scratch
1111int ClpPredictorCorrector::createSolution()
1112{
1113  int numberTotal = numberRows_+numberColumns_;
1114  int iColumn;
1115  double tolerance = primalTolerance();
1116  for (iColumn=0;iColumn<numberTotal;iColumn++) {
1117    if (upper_[iColumn]-lower_[iColumn]>tolerance) 
1118      clearFixed(iColumn);
1119    else
1120      setFixed(iColumn);
1121  }
1122  double initialValue =1.0e-12;
1123
1124  double maximumObjective=1.0e-12;
1125  double objectiveNorm2=0.0;
1126  getNorms(cost_,numberTotal,maximumObjective,objectiveNorm2);
1127  if (maximumObjective<1.0e-12) {
1128    maximumObjective=1.0e-12;
1129  } 
1130  objectiveNorm2=sqrt(objectiveNorm2)/(double) numberTotal;
1131  objectiveNorm_=maximumObjective;
1132  scaleFactor_=1.0;
1133  if (maximumObjective>0.0) {
1134    if (maximumObjective<1.0) {
1135      scaleFactor_=maximumObjective;
1136    } else if (maximumObjective>1.0e4) {
1137      scaleFactor_=maximumObjective/1.0e4;
1138    } 
1139  } 
1140  if (scaleFactor_!=1.0) {
1141    objectiveNorm2*=scaleFactor_;
1142    multiplyAdd(NULL,numberTotal,0.0,cost_,1.0/scaleFactor_);
1143    objectiveNorm_=maximumObjective/scaleFactor_;
1144  } 
1145  baseObjectiveNorm_=objectiveNorm_;
1146  //accumulate fixed in dj region (as spare)
1147  //accumulate primal solution in primal region
1148  //DZ in lowerDual
1149  //DW in upperDual
1150  double infiniteCheck=1.0e40;
1151  //double     fakeCheck=1.0e10;
1152  //use deltaX region for work region
1153  for (iColumn=0;iColumn<numberTotal;iColumn++) {
1154    double primalValue = solution_[iColumn];
1155    clearFlagged(iColumn);
1156    clearFixedOrFree(iColumn);
1157    clearLowerBound(iColumn);
1158    clearUpperBound(iColumn);
1159    clearFakeLower(iColumn);
1160    clearFakeUpper(iColumn);
1161    if (!fixed(iColumn)) {
1162      dj_[iColumn]=0.0;
1163      diagonal_[iColumn]=1.0;
1164      deltaX_[iColumn]=1.0;
1165      double lowerValue=lower_[iColumn];
1166      double upperValue=upper_[iColumn];
1167#if 0
1168      //fake it
1169      if (lowerValue<-fakeCheck&&upperValue>fakeCheck) {
1170        lowerValue=-1.0e5;
1171        upperValue=1.0e5;
1172        lower_[iColumn]=lowerValue;
1173        upper_[iColumn]=upperValue;
1174        std::cout<<"faking free variable "<<iColumn<<std::endl;
1175      }
1176#endif
1177      if (lowerValue>-infiniteCheck) {
1178        if (upperValue<infiniteCheck) {
1179          //upper and lower bounds
1180          setLowerBound(iColumn);
1181          setUpperBound(iColumn);
1182          if (lowerValue>=0.0) {
1183            solution_[iColumn]=lowerValue;
1184          } else if (upperValue<=0.0) {
1185            solution_[iColumn]=upperValue;
1186          } else {
1187            solution_[iColumn]=0.0;
1188          } 
1189        } else {
1190          //just lower bound
1191          setLowerBound(iColumn);
1192          if (lowerValue>=0.0) {
1193            solution_[iColumn]=lowerValue;
1194          } else {
1195            solution_[iColumn]=0.0;
1196          } 
1197        } 
1198      } else {
1199        if (upperValue<infiniteCheck) {
1200          //just upper bound
1201          setUpperBound(iColumn);
1202          if (upperValue<=0.0) {
1203            solution_[iColumn]=upperValue;
1204          } else {
1205            solution_[iColumn]=0.0;
1206          } 
1207        } else {
1208          //free
1209          setFixedOrFree(iColumn);
1210          solution_[iColumn]=0.0;
1211          //std::cout<<" free "<<i<<std::endl;
1212        } 
1213      } 
1214    } else {
1215      setFlagged(iColumn);
1216      setFixedOrFree(iColumn);
1217      setLowerBound(iColumn);
1218      setUpperBound(iColumn);
1219      dj_[iColumn]=primalValue;;
1220      solution_[iColumn]=lower_[iColumn];
1221      diagonal_[iColumn]=0.0;
1222      deltaX_[iColumn]=0.0;
1223    } 
1224  } 
1225  //   modify fixed RHS
1226  multiplyAdd(dj_+numberColumns_,numberRows_,1.0,rhsFixRegion_,0.0);
1227  //   create plausible RHS?
1228  matrix_->times(-1.0,dj_,rhsFixRegion_);
1229  multiplyAdd(solution_+numberColumns_,numberRows_,-1.0,errorRegion_,0.0);
1230  matrix_->times(1.0,solution_,errorRegion_);
1231  rhsNorm_=maximumAbsElement(errorRegion_,numberRows_);
1232  if (rhsNorm_<1.0) {
1233    rhsNorm_=1.0;
1234  } 
1235  int * rowsDropped = new int [numberRows_];
1236  int returnCode=cholesky_->factorize(diagonal_,rowsDropped);
1237  if (returnCode==-1) {
1238    printf("Out of memory\n");
1239    problemStatus_=4;
1240    return -1;
1241  }
1242  if (cholesky_->status()) {
1243    std::cout << "singular on initial cholesky?" <<std::endl;
1244    cholesky_->resetRowsDropped();
1245    //cholesky_->factorize(rowDropped_);
1246    //if (cholesky_->status()) {
1247      //std::cout << "bad cholesky??? (after retry)" <<std::endl;
1248      //abort();
1249    //}
1250  } 
1251  delete [] rowsDropped;
1252#ifndef KKT
1253  cholesky_->solve(errorRegion_);
1254  //create information for solution
1255  multiplyAdd(errorRegion_,numberRows_,-1.0,deltaX_+numberColumns_,0.0);
1256  CoinZeroN(deltaX_,numberColumns_);
1257  matrix_->transposeTimes(1.0,errorRegion_,deltaX_);
1258#else
1259  // reverse sign on solution
1260  multiplyAdd(NULL,numberRows_+numberColumns_,0.0,solution_,-1.0);
1261  solveSystem(deltaX_,errorRegion_,solution_,NULL,NULL,NULL,false);
1262#endif
1263  //do reduced costs
1264  CoinZeroN(dj_+numberColumns_,numberRows_);
1265  for ( iColumn=0;iColumn<numberColumns_;iColumn++) {
1266    double value=cost_[iColumn];
1267    if (lowerBound(iColumn)) {
1268      value+=linearPerturbation_;
1269    } 
1270    if (upperBound(iColumn)) {
1271      value-=linearPerturbation_;
1272    } 
1273    dj_[iColumn]=value;
1274  } 
1275  initialValue=1.0e1;
1276  if (rhsNorm_*1.0e-2>initialValue) {
1277    initialValue=rhsNorm_*1.0e-2;
1278  } 
1279  double smallestBoundDifference=COIN_DBL_MAX;
1280  double * fakeSolution = deltaX_;
1281  for ( iColumn=0;iColumn<numberTotal;iColumn++) {
1282    if (!flagged(iColumn)) {
1283      if (lower_[iColumn]-fakeSolution[iColumn]>initialValue) {
1284        initialValue=lower_[iColumn]-fakeSolution[iColumn];
1285      } 
1286      if (fakeSolution[iColumn]-upper_[iColumn]>initialValue) {
1287        initialValue=fakeSolution[iColumn]-upper_[iColumn];
1288      } 
1289      if (upper_[iColumn]-lower_[iColumn]<smallestBoundDifference) {
1290        smallestBoundDifference=upper_[iColumn]-lower_[iColumn];
1291      } 
1292    } 
1293  } 
1294  solutionNorm_=1.0e-12;
1295  handler_->message(CLP_BARRIER_SAFE,messages_)
1296    <<initialValue<<objectiveNorm_
1297    <<CoinMessageEol;
1298  int strategy=0;
1299  double extra=1.0e-10;
1300  double largeGap=1.0e15;
1301  double safeObjectiveValue=2.0*objectiveNorm_;
1302  double safeFree=1.0e-1*initialValue;
1303  //safeFree=1.0;
1304  //printf("normal safe dual value of %g, primal value of %g\n",
1305  // safeObjectiveValue,initialValue);
1306  //safeObjectiveValue=max(2.0,1.0e-1*safeObjectiveValue);
1307  //initialValue=max(100.0,1.0e-1*initialValue);;
1308  //printf("temp safe dual value of %g, primal value of %g\n",
1309  // safeObjectiveValue,initialValue);
1310  double zwLarge=1.0e2*initialValue;
1311  //zwLarge=1.0e40;
1312  for ( iColumn=0;iColumn<numberTotal;iColumn++) {
1313    if (!flagged(iColumn)) {
1314      double lowerValue=lower_[iColumn];
1315      double upperValue=upper_[iColumn];
1316      double objectiveValue = cost_[iColumn];
1317      double newValue;
1318      double low;
1319      double high;
1320      double randomZ =0.1*CoinDrand48()+0.9;
1321      double randomW =0.1*CoinDrand48()+0.9;
1322      double setPrimal=initialValue;
1323      if (strategy==0) {
1324        randomZ=1.0;
1325        randomW=1.0;
1326      } 
1327      if (lowerBound(iColumn)) {
1328        if (upperBound(iColumn)) {
1329          //upper and lower bounds
1330          if (upperValue-lowerValue>2.0*setPrimal) {
1331            double fakeValue=fakeSolution[iColumn];
1332            if (fakeValue<lowerValue+setPrimal) {
1333              fakeValue=lowerValue+setPrimal;
1334            } 
1335            if (fakeValue>upperValue-setPrimal) {
1336              fakeValue=upperValue-setPrimal;
1337            } 
1338            newValue=fakeValue;
1339            low=fakeValue-lowerValue;
1340            high=upperValue-fakeValue;
1341          } else {
1342            newValue=0.5*(upperValue+lowerValue);
1343            low=setPrimal;
1344            high=setPrimal;
1345          } 
1346          low *= randomZ;
1347          high *= randomW;
1348          double s = low+extra;
1349          double ratioZ;
1350          if (s<zwLarge) {
1351            ratioZ=1.0;
1352          } else {
1353            ratioZ=sqrt(zwLarge/s);
1354          } 
1355          double t = high+extra;
1356          double ratioW;
1357          if (t<zwLarge) {
1358            ratioW=1.0;
1359          } else {
1360            ratioW=sqrt(zwLarge/t);
1361          } 
1362          //modify s and t
1363          if (s>largeGap) {
1364            s=largeGap;
1365          } 
1366          if (t>largeGap) {
1367            t=largeGap;
1368          } 
1369          //modify if long long way away from bound
1370          if (objectiveValue>=0.0) {
1371            zVec_[iColumn]=objectiveValue + safeObjectiveValue*ratioZ;
1372            tVec_[iColumn]=safeObjectiveValue*ratioW;
1373          } else {
1374            zVec_[iColumn]=safeObjectiveValue*ratioZ;
1375            tVec_[iColumn]=-objectiveValue + safeObjectiveValue*ratioW;
1376          } 
1377          diagonal_[iColumn] = (t*s)/(s*tVec_[iColumn]+t*zVec_[iColumn]);
1378        } else {
1379          //just lower bound
1380          double fakeValue=fakeSolution[iColumn];
1381          if (fakeValue<lowerValue+setPrimal) {
1382            fakeValue=lowerValue+setPrimal;
1383          } 
1384          newValue=fakeValue;
1385          low=fakeValue-lowerValue;
1386          low *= randomZ;
1387          high=0.0;
1388          double s = low+extra;
1389          double ratioZ;
1390          if (s<zwLarge) {
1391            ratioZ=1.0;
1392          } else {
1393            ratioZ=sqrt(zwLarge/s);
1394          } 
1395          //modify s
1396          if (s>largeGap) {
1397            s=largeGap;
1398          } 
1399          if (objectiveValue>=0.0) {
1400            zVec_[iColumn]=objectiveValue + safeObjectiveValue*ratioZ;
1401            tVec_[iColumn]=0.0;
1402          } else {
1403            zVec_[iColumn]=safeObjectiveValue*ratioZ;
1404            tVec_[iColumn]=0.0;
1405          } 
1406          diagonal_[iColumn] = s/zVec_[iColumn];
1407        } 
1408      } else {
1409        if (upperBound(iColumn)) {
1410          //just upper bound
1411          double fakeValue=fakeSolution[iColumn];
1412          if (fakeValue>upperValue-setPrimal) {
1413            fakeValue=upperValue-setPrimal;
1414          } 
1415          newValue=fakeValue;
1416          low=0.0;
1417          high=upperValue-fakeValue;
1418          high*=randomW;
1419          double t = high+extra;
1420          double ratioW;
1421          if (t<zwLarge) {
1422            ratioW=1.0;
1423          } else {
1424            ratioW=sqrt(zwLarge/t);
1425          } 
1426          //modify t
1427          if (t>largeGap) {
1428            t=largeGap;
1429          } 
1430          if (objectiveValue>=0.0) {
1431            zVec_[iColumn]=0.0;
1432            tVec_[iColumn]=safeObjectiveValue*ratioW;
1433          } else {
1434            zVec_[iColumn]=0.0;
1435            tVec_[iColumn]=-objectiveValue + safeObjectiveValue*ratioW;
1436          } 
1437          diagonal_[iColumn] =  t/tVec_[iColumn];
1438        } else {
1439          //free
1440          zVec_[iColumn]=safeObjectiveValue;
1441          tVec_[iColumn]=safeObjectiveValue;
1442          newValue=fakeSolution[iColumn];
1443          if (newValue>=0.0) {
1444            if (newValue<safeFree) {
1445              newValue=safeFree;
1446            } 
1447          } else {
1448            if (newValue>-safeFree) {
1449              newValue=-safeFree;
1450            } 
1451          } 
1452          if (fabs(newValue)>1.0) {
1453            diagonal_[iColumn]=fabs(newValue)/(zVec_[iColumn]+tVec_[iColumn]);
1454          } else {
1455            diagonal_[iColumn]=1.0/(zVec_[iColumn]+tVec_[iColumn]);
1456          } 
1457          low=0.0;
1458          high=0.0;
1459        } 
1460      } 
1461      lowerSlack_[iColumn]=low;
1462      upperSlack_[iColumn]=high;
1463      solution_[iColumn]=newValue;
1464    } else {
1465      lowerSlack_[iColumn]=0.0;
1466      upperSlack_[iColumn]=0.0;
1467      solution_[iColumn]=lower_[iColumn];
1468      zVec_[iColumn]=0.0;
1469      tVec_[iColumn]=0.0;
1470      diagonal_[iColumn]=0.0;
1471    } 
1472  } 
1473  solutionNorm_ =  maximumAbsElement(solution_,numberTotal);
1474  return 0;
1475}
1476// complementarityGap.  Computes gap
1477//phase 0=as is , 1 = after predictor , 2 after corrector
1478double ClpPredictorCorrector::complementarityGap(int & numberComplementarityPairs,
1479                                                 int & numberComplementarityItems,
1480                                                 const int phase)
1481{
1482  double gap=0.0;
1483  //seems to be same coding for phase = 1 or 2
1484  numberComplementarityPairs=0;
1485  numberComplementarityItems=0;
1486  double toleranceGap=0.0;
1487  double largestGap=0.0;
1488  double smallestGap=COIN_DBL_MAX;
1489  //seems to be same coding for phase = 1 or 2
1490  int numberNegativeGaps=0;
1491  double sumNegativeGap=0.0;
1492  double largeGap=1.0e2*solutionNorm_;
1493  if (largeGap<1.0e10) {
1494    largeGap=1.0e10;
1495  } 
1496  largeGap=1.0e30;
1497  double dualTolerance =  dblParam_[ClpDualTolerance];
1498  double primalTolerance =  dblParam_[ClpPrimalTolerance];
1499  dualTolerance=dualTolerance/scaleFactor_;
1500  double * zVec = zVec_;
1501  double * tVec = tVec_;
1502  double * primal = solution_;
1503  double * lower = lower_;
1504  double * upper = upper_;
1505  double * lowerSlack = lowerSlack_;
1506  double * upperSlack = upperSlack_;
1507  double * deltaZ = deltaZ_;
1508  double * deltaT = deltaT_;
1509  double * deltaX = deltaX_;
1510  for (int iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
1511    if (!fixedOrFree(iColumn)) {
1512      numberComplementarityPairs++;
1513      //can collapse as if no lower bound both zVec and deltaZ 0.0
1514      if (lowerBound(iColumn)) {
1515        numberComplementarityItems++;
1516        double dualValue;
1517        double primalValue;
1518        if (!phase) {
1519          dualValue=zVec[iColumn];
1520          primalValue=lowerSlack[iColumn];
1521        } else {
1522          double change;
1523          if (!fakeLower(iColumn)) {
1524            change =primal[iColumn]+deltaX[iColumn]-lowerSlack[iColumn]-lower[iColumn];
1525          } else {
1526            change =deltaX[iColumn];
1527          } 
1528          dualValue=zVec[iColumn]+actualDualStep_*deltaZ[iColumn];
1529          primalValue=lowerSlack[iColumn]+actualPrimalStep_*change;
1530        } 
1531        //reduce primalValue
1532        if (primalValue>largeGap) {
1533          primalValue=largeGap;
1534        } 
1535        double gapProduct=dualValue*primalValue;
1536        if (gapProduct<0.0) {
1537          //cout<<"negative gap component "<<iColumn<<" "<<dualValue<<" "<<
1538              //primalValue<<endl;
1539          numberNegativeGaps++;
1540          sumNegativeGap-=gapProduct;
1541          gapProduct=0.0;
1542        } 
1543        gap+=gapProduct;
1544        if (gapProduct>largestGap) {
1545          largestGap=gapProduct;
1546        }
1547        smallestGap = min(smallestGap,gapProduct);
1548        if (dualValue>dualTolerance&&primalValue>primalTolerance) {
1549          toleranceGap+=dualValue*primalValue;
1550        } 
1551      } 
1552      if (upperBound(iColumn)) {
1553        numberComplementarityItems++;
1554        double dualValue;
1555        double primalValue;
1556        if (!phase) {
1557          dualValue=tVec[iColumn];
1558          primalValue=upperSlack[iColumn];
1559        } else {
1560          double change;
1561          if (!fakeUpper(iColumn)) {
1562            change =upper[iColumn]-primal[iColumn]-deltaX[iColumn]-upperSlack[iColumn];
1563          } else {
1564            change =-deltaX[iColumn];
1565          } 
1566          dualValue=tVec[iColumn]+actualDualStep_*deltaT[iColumn];
1567          primalValue=upperSlack[iColumn]+actualPrimalStep_*change;
1568        } 
1569        //reduce primalValue
1570        if (primalValue>largeGap) {
1571          primalValue=largeGap;
1572        } 
1573        double gapProduct=dualValue*primalValue;
1574        if (gapProduct<0.0) {
1575          //cout<<"negative gap component "<<iColumn<<" "<<dualValue<<" "<<
1576              //primalValue<<endl;
1577          numberNegativeGaps++;
1578          sumNegativeGap-=gapProduct;
1579          gapProduct=0.0;
1580        } 
1581        gap+=gapProduct;
1582        if (gapProduct>largestGap) {
1583          largestGap=gapProduct;
1584        } 
1585        if (dualValue>dualTolerance&&primalValue>primalTolerance) {
1586          toleranceGap+=dualValue*primalValue;
1587        } 
1588      } 
1589    } 
1590  } 
1591  if (!phase&&numberNegativeGaps) {
1592      handler_->message(CLP_BARRIER_NEGATIVE_GAPS,messages_)
1593    <<numberNegativeGaps<<sumNegativeGap
1594    <<CoinMessageEol;
1595  } 
1596 
1597  //in case all free!
1598  if (!numberComplementarityPairs) {
1599    numberComplementarityPairs=1;
1600  } 
1601  //printf("gap %g - smallest %g, largest %g, pairs %d\n",
1602  // gap,smallestGap,largestGap,numberComplementarityPairs);
1603  return gap;
1604}
1605// setupForSolve.
1606//phase 0=affine , 1 = corrector , 2 = primal-dual
1607void ClpPredictorCorrector::setupForSolve(const int phase)
1608{
1609  double extra =eExtra;
1610  double * zVec = zVec_;
1611  double * tVec = tVec_;
1612  double * primal = solution_;
1613  double * dj = dj_;
1614  double * lower = lower_;
1615  double * upper = upper_;
1616  double * lowerSlack = lowerSlack_;
1617  double * upperSlack = upperSlack_;
1618  double * workArray = workArray_;
1619
1620  int iColumn;
1621#ifdef SOME_DEBUG
1622  printf("phase %d in setupForSolve, mu %g\n",phase,mu_);
1623#endif
1624  switch (phase) {
1625  case 0:
1626    CoinMemcpyN(errorRegion_,numberRows_,rhsB_);
1627    for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
1628      rhsC_[iColumn]=0.0;
1629      rhsU_[iColumn]=0.0;
1630      rhsL_[iColumn]=0.0;
1631      rhsZ_[iColumn]=0.0;
1632      rhsT_[iColumn]=0.0;
1633      if (!flagged(iColumn)) {
1634        rhsC_[iColumn] = dj[iColumn]-zVec[iColumn]+tVec[iColumn];
1635        if (lowerBound(iColumn)) {
1636          if (!fakeLower(iColumn)) {
1637            rhsZ_[iColumn] = -zVec[iColumn]*(lowerSlack[iColumn]+extra);
1638            rhsL_[iColumn] = primal[iColumn]-lowerSlack[iColumn]-lower[iColumn];
1639          } 
1640        } 
1641        if (upperBound(iColumn)) {
1642          if (!fakeUpper(iColumn)) {
1643            rhsT_[iColumn] = -tVec[iColumn]*(upperSlack[iColumn]+extra);
1644            rhsU_[iColumn] = -(primal[iColumn] + upperSlack[iColumn] - upper[iColumn]); 
1645          } 
1646        }
1647      }
1648    } 
1649    break;
1650  case 1:
1651    // could be stored in delta2?
1652    for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
1653      rhsC_[iColumn]=0.0;
1654      //rhsU_[iColumn]=0.0;
1655      //rhsL_[iColumn]=0.0;
1656      rhsZ_[iColumn]=0.0;
1657      rhsT_[iColumn]=0.0;
1658      if (!flagged(iColumn)) {
1659        rhsC_[iColumn] = dj[iColumn]-zVec[iColumn]+tVec[iColumn];
1660        if (lowerBound(iColumn)) {
1661          rhsZ_[iColumn] = mu_ -zVec[iColumn]*(lowerSlack[iColumn]+extra)
1662            - deltaZ_[iColumn]*deltaX_[iColumn];
1663          // To bring in line with OSL
1664          //rhsZ_[iColumn] -= deltaZ_[iColumn]*rhsL_[iColumn];
1665        } 
1666        if (upperBound(iColumn)) {
1667          rhsT_[iColumn] = mu_ -tVec[iColumn]*(upperSlack[iColumn]+extra)
1668            +deltaT_[iColumn]*deltaX_[iColumn];
1669          // To bring in line with OSL
1670          //rhsT_[iColumn] += deltaT_[iColumn]*rhsU_[iColumn];
1671        } 
1672      } 
1673    } 
1674    break;
1675  case 2:
1676    CoinMemcpyN(errorRegion_,numberRows_,rhsB_);
1677    for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
1678      rhsC_[iColumn]=0.0;
1679      rhsU_[iColumn]=0.0;
1680      rhsL_[iColumn]=0.0;
1681      rhsZ_[iColumn]=0.0;
1682      rhsT_[iColumn]=0.0;
1683      if (!flagged(iColumn)) {
1684        rhsC_[iColumn] = dj[iColumn]-zVec[iColumn]+tVec[iColumn];
1685        if (lowerBound(iColumn)) {
1686          if (!fakeLower(iColumn)) {
1687            rhsZ_[iColumn] = mu_ - zVec[iColumn]*(lowerSlack[iColumn]+extra);
1688            rhsL_[iColumn] = primal[iColumn]-lowerSlack[iColumn]-lower[iColumn];
1689          } 
1690        } 
1691        if (upperBound(iColumn)) {
1692          if (!fakeUpper(iColumn)) {
1693            rhsT_[iColumn] = mu_ - tVec[iColumn]*(upperSlack[iColumn]+extra);
1694            rhsU_[iColumn] = -(primal[iColumn] + upperSlack[iColumn] - upper[iColumn]); 
1695          } 
1696        }
1697      }
1698    } 
1699    break;
1700  case 3:
1701    {
1702      double minBeta = 0.1*mu_;
1703      double maxBeta = 10.0*mu_;
1704      double dualStep = min(1.0,actualDualStep_+0.1);
1705      double primalStep = min(1.0,actualPrimalStep_+0.1);
1706#ifdef SOME_DEBUG
1707      printf("good complementarity range %g to %g\n",minBeta,maxBeta);
1708#endif
1709      //minBeta=0.0;
1710      //maxBeta=COIN_DBL_MAX;
1711      for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
1712        if (!flagged(iColumn)) {
1713          if (lowerBound(iColumn)) {
1714            double change = rhsL_[iColumn] + deltaX_[iColumn];
1715            double dualValue=zVec[iColumn]+dualStep*deltaZ_[iColumn];
1716            double primalValue=lowerSlack[iColumn]+primalStep*change;
1717            double gapProduct=dualValue*primalValue;
1718            if (gapProduct>0.0&&dualValue<0.0)
1719              gapProduct = - gapProduct;
1720#ifdef FULL_DEBUG
1721            delta2Z_[iColumn]=gapProduct;
1722            if (delta2Z_[iColumn]<minBeta||delta2Z_[iColumn]>maxBeta)
1723              printf("lower %d primal %g, dual %g, gap %g\n",
1724                     iColumn,primalValue,dualValue,gapProduct);
1725#endif
1726            double value= 0.0;
1727            if (gapProduct<minBeta) {
1728              value= 2.0*(minBeta-gapProduct);
1729              value= (mu_-gapProduct);
1730              value= (minBeta-gapProduct);
1731              assert (value>0.0);
1732            } else if (gapProduct>maxBeta) {
1733              value= max(maxBeta-gapProduct,-maxBeta);
1734              assert (value<0.0);
1735            }
1736            //#define AGAIN
1737#ifdef AGAIN
1738            //rhsZ_[iColumn] = mu_ -zVec[iColumn]*(lowerSlack[iColumn]+extra)
1739            //- deltaZ_[iColumn]*deltaX_[iColumn];
1740            // To bring in line with OSL
1741            rhsZ_[iColumn] += deltaZ_[iColumn]*rhsL_[iColumn];
1742#endif
1743            rhsZ_[iColumn] += value;
1744          } 
1745          if (upperBound(iColumn)) {
1746            double change = rhsU_[iColumn]-deltaX_[iColumn];
1747            double dualValue=tVec[iColumn]+dualStep*deltaT_[iColumn];
1748            double primalValue=upperSlack[iColumn]+primalStep*change;
1749            double gapProduct=dualValue*primalValue;
1750            if (gapProduct>0.0&&dualValue<0.0)
1751              gapProduct = - gapProduct;
1752#ifdef FULL_DEBUG
1753            delta2T_[iColumn]=gapProduct;
1754            if (delta2T_[iColumn]<minBeta||delta2T_[iColumn]>maxBeta)
1755              printf("upper %d primal %g, dual %g, gap %g\n",
1756                     iColumn,primalValue,dualValue,gapProduct);
1757#endif
1758            double value= 0.0;
1759            if (gapProduct<minBeta) {
1760              value= (minBeta-gapProduct);
1761              assert (value>0.0);
1762            } else if (gapProduct>maxBeta) {
1763              value= max(maxBeta-gapProduct,-maxBeta);
1764              assert (value<0.0);
1765            }
1766#ifdef AGAIN
1767            //rhsT_[iColumn] = mu_ -tVec[iColumn]*(upperSlack[iColumn]+extra)
1768            //+deltaT_[iColumn]*deltaX_[iColumn];
1769            // To bring in line with OSL
1770            rhsT_[iColumn] -= deltaT_[iColumn]*rhsU_[iColumn];
1771#endif
1772            rhsT_[iColumn] += value;
1773          } 
1774        } 
1775      }
1776    }
1777    break;
1778  } /* endswitch */
1779#ifndef KKT
1780  for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
1781    double value = rhsC_[iColumn];
1782    if (lowerBound(iColumn)) {
1783      value -= rhsZ_[iColumn]/(lowerSlack[iColumn]+extra);
1784      value += zVec[iColumn]*rhsL_[iColumn]/(lowerSlack[iColumn]+extra);
1785    }
1786    if (upperBound(iColumn)) {
1787      value += rhsT_[iColumn]/(upperSlack[iColumn]+extra);
1788      value -= tVec[iColumn]*rhsU_[iColumn]/(upperSlack[iColumn]+extra);
1789    }
1790    workArray[iColumn]=diagonal_[iColumn]*value;
1791  } 
1792#else
1793  for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
1794    double value = rhsC_[iColumn];
1795    if (lowerBound(iColumn)) {
1796      value -= rhsZ_[iColumn]/(lowerSlack[iColumn]+extra);
1797      value += zVec[iColumn]*rhsL_[iColumn]/(lowerSlack[iColumn]+extra);
1798    }
1799    if (upperBound(iColumn)) {
1800      value += rhsT_[iColumn]/(upperSlack[iColumn]+extra);
1801      value -= tVec[iColumn]*rhsU_[iColumn]/(upperSlack[iColumn]+extra);
1802    }
1803      workArray[iColumn]=value;
1804  }
1805#endif
1806}
1807//method: sees if looks plausible change in complementarity
1808bool ClpPredictorCorrector::checkGoodMove(const bool doCorrector,double & bestNextGap)
1809{
1810  const double beta3 = 0.99997;
1811  bool goodMove=false;
1812  int nextNumber;
1813  int nextNumberItems;
1814  int numberTotal = numberRows_+numberColumns_;
1815  double returnGap=bestNextGap;
1816  double nextGap=complementarityGap(nextNumber,nextNumberItems,2);
1817  if (nextGap>bestNextGap&&nextGap>-0.99*complementarityGap_) {
1818#ifdef SOME_DEBUG
1819    printf("checkGood phase 1 next gap %g, phase 0 %g, old gap %g\n",
1820           nextGap,bestNextGap,complementarityGap_);
1821#endif
1822    return false;
1823  } else {
1824    returnGap=nextGap;
1825  }
1826  double step;
1827  if (actualDualStep_>actualPrimalStep_) {
1828    step=actualDualStep_;
1829  } else {
1830    step=actualPrimalStep_;
1831  } 
1832  double testValue=1.0-step*(1.0-beta3);
1833  //testValue=0.0;
1834  testValue*=complementarityGap_;
1835  if (nextGap<testValue) {
1836    //std::cout <<"predicted duality gap "<<nextGap<<std::endl;
1837    goodMove=true;
1838  } else if(doCorrector) {
1839    //if (actualDualStep_<actualPrimalStep_) {
1840      //step=actualDualStep_;
1841    //} else {
1842      //step=actualPrimalStep_;
1843    //}
1844    double gap = bestNextGap;
1845    goodMove=checkGoodMove2(step,gap);
1846    if (goodMove)
1847      returnGap=gap;
1848  } else {
1849    goodMove=true;
1850  } 
1851  if (!goodMove) {
1852    //try smaller of two
1853    if (actualDualStep_<actualPrimalStep_) {
1854      step=actualDualStep_;
1855    } else {
1856      step=actualPrimalStep_;
1857    } 
1858    if (step>1.0) {
1859      step=1.0;
1860    } 
1861    actualPrimalStep_=step;
1862    actualDualStep_=step;
1863    while (!goodMove) {
1864      double gap = bestNextGap;
1865      goodMove=checkGoodMove2(step,gap);
1866      if (goodMove)
1867        returnGap=gap;
1868      if (step<1.0e-10) {
1869        break;
1870      } 
1871      step*=0.5;
1872      actualPrimalStep_=step;
1873      actualDualStep_=step;
1874    } /* endwhile */
1875    if (doCorrector) {
1876      //say bad move if both small
1877      if (numberIterations_&1) {
1878        if (actualPrimalStep_<1.0e-2&&actualDualStep_<1.0e-2) {
1879          goodMove=false;
1880        } 
1881      } else {
1882        if (actualPrimalStep_<1.0e-5&&actualDualStep_<1.0e-5) {
1883          goodMove=false;
1884        } 
1885        if (actualPrimalStep_*actualDualStep_<1.0e-20) {
1886          goodMove=false;
1887        } 
1888      } 
1889    } 
1890  } 
1891  if (goodMove) {
1892    //compute delta in objectives
1893    double deltaObjectivePrimal=0.0;
1894    double deltaObjectiveDual=
1895      innerProduct(deltaY_,numberRows_,
1896                   rhsFixRegion_);
1897    double error=0.0;
1898    double * workArray = workArray_;
1899    CoinZeroN(workArray,numberColumns_);
1900    CoinMemcpyN(deltaY_,numberRows_,workArray+numberColumns_);
1901    matrix_->transposeTimes(-1.0,deltaY_,workArray);
1902    double * deltaZ = deltaZ_;
1903    double * deltaT = deltaT_;
1904    double * lower = lower_;
1905    double * upper = upper_;
1906    //direction vector in deltaX
1907    double * deltaX = deltaX_;
1908    double * cost = cost_;
1909    //double sumPerturbCost=0.0;
1910    for (int iColumn=0;iColumn<numberTotal;iColumn++) {
1911      if (!flagged(iColumn)) {
1912        if (lowerBound(iColumn)) {
1913          //sumPerturbCost+=deltaX[iColumn];
1914          deltaObjectiveDual+=deltaZ[iColumn]*lower[iColumn];
1915        } 
1916        if (upperBound(iColumn)) {
1917          //sumPerturbCost-=deltaX[iColumn];
1918          deltaObjectiveDual-=deltaT[iColumn]*upper[iColumn];
1919        } 
1920        double change = fabs(workArray[iColumn]-deltaZ[iColumn]+deltaT[iColumn]);
1921        error = max (change,error);
1922      } 
1923      deltaObjectivePrimal += cost[iColumn] * deltaX[iColumn];
1924    } 
1925    //deltaObjectivePrimal+=sumPerturbCost*linearPerturbation_;
1926    double testValue;
1927    if (error>0.0) {
1928      testValue=1.0e1*maximumDualError_/error;
1929    } else {
1930      testValue=1.0e1;
1931    } 
1932    if (testValue<actualDualStep_) {
1933      handler_->message(CLP_BARRIER_REDUCING,messages_)
1934      <<"dual"<<actualDualStep_
1935      << testValue
1936      <<CoinMessageEol;
1937      actualDualStep_=testValue;
1938    } 
1939  } 
1940  if (maximumRHSError_<1.0e1*solutionNorm_*primalTolerance()
1941                            &&maximumRHSChange_>1.0e-16*solutionNorm_) {
1942    //check change in AX not too much
1943    //??? could be dropped row going infeasible
1944    double ratio = 1.0e1*maximumRHSError_/maximumRHSChange_;
1945    if (ratio<actualPrimalStep_) {
1946      handler_->message(CLP_BARRIER_REDUCING,messages_)
1947      <<"primal"<<actualPrimalStep_
1948      <<ratio
1949      <<CoinMessageEol;
1950      if (ratio>1.0e-6) {
1951        actualPrimalStep_=ratio;
1952      } else {
1953        actualPrimalStep_=ratio;
1954        //std::cout <<"sign we should be stopping"<<std::endl;
1955      } 
1956    } 
1957  }
1958  if (goodMove)
1959    bestNextGap=returnGap;
1960  return goodMove;
1961}
1962//:  checks for one step size
1963bool ClpPredictorCorrector::checkGoodMove2(const double move,double & bestNextGap)
1964{
1965  double complementarityMultiplier =1.0/numberComplementarityPairs_;
1966  const double gamma = 1.0e-8;
1967  const double gammap = 1.0e-8;
1968  const double gammad = 1.0e-8;
1969  int nextNumber;
1970  int nextNumberItems;
1971  double nextGap=complementarityGap(nextNumber,nextNumberItems,2);
1972  if (nextGap>bestNextGap)
1973    return false;
1974  double lowerBoundGap = gamma*nextGap*complementarityMultiplier;
1975  bool goodMove=true;
1976  double * deltaZ = deltaZ_;
1977  double * deltaT = deltaT_;
1978  double * deltaSL = deltaSL_;
1979  double * deltaSU = deltaSU_;
1980  double * zVec = zVec_;
1981  double * tVec = tVec_;
1982  double * lowerSlack = lowerSlack_;
1983  double * upperSlack = upperSlack_;
1984  for (int iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
1985    if (!flagged(iColumn)) {
1986      if (lowerBound(iColumn)) {
1987        double part1=lowerSlack[iColumn]+actualPrimalStep_*deltaSL[iColumn];
1988        double part2=zVec[iColumn]+actualDualStep_*deltaZ[iColumn];
1989        if (part1*part2<lowerBoundGap) {
1990          goodMove=false;
1991          break;
1992        } 
1993      } 
1994      if (upperBound(iColumn)) {
1995        double part1=upperSlack[iColumn]+actualPrimalStep_*deltaT[iColumn];
1996        double part2=tVec[iColumn]+actualDualStep_*deltaSU[iColumn];
1997        if (part1*part2<lowerBoundGap) {
1998          goodMove=false;
1999          break;
2000        } 
2001      } 
2002    } 
2003  } 
2004  //      Satisfy g_p(alpha)?
2005  if (rhsNorm_>solutionNorm_) {
2006    solutionNorm_=rhsNorm_;
2007  } 
2008  double errorCheck=maximumRHSError_/solutionNorm_;
2009  if (errorCheck<maximumBoundInfeasibility_) {
2010    errorCheck=maximumBoundInfeasibility_;
2011  } 
2012  //scale
2013  if ((1.0-move)*errorCheck>primalTolerance()) {
2014    if (nextGap<gammap*(1.0-move)*errorCheck) {
2015      goodMove=false;
2016    } 
2017  } 
2018  //      Satisfy g_d(alpha)?
2019  errorCheck=maximumDualError_/objectiveNorm_;
2020  if ((1.0-move)*errorCheck>dualTolerance()) {
2021    if (nextGap<gammad*(1.0-move)*errorCheck) {
2022      goodMove=false;
2023    } 
2024  }
2025  if (goodMove)
2026    bestNextGap=nextGap;
2027  return goodMove;
2028}
2029// updateSolution.  Updates solution at end of iteration
2030//returns number fixed
2031int ClpPredictorCorrector::updateSolution()
2032{
2033  //update pi
2034  multiplyAdd(deltaY_,numberRows_,actualDualStep_,dual_,1.0);
2035  CoinZeroN(errorRegion_,numberRows_);
2036  CoinZeroN(rhsFixRegion_,numberRows_);
2037  double maximumRhsInfeasibility=0.0;
2038  double maximumBoundInfeasibility=0.0;
2039  double maximumDualError=1.0e-12;
2040  double primalObjectiveValue=0.0;
2041  double dualObjectiveValue=0.0;
2042  double solutionNorm=1.0e-12;
2043  int numberKilled=0;
2044  double freeMultiplier=1.0e6;
2045  double trueNorm =diagonalNorm_/diagonalScaleFactor_;
2046  if (freeMultiplier<trueNorm) {
2047    freeMultiplier=trueNorm;
2048  } 
2049  if (freeMultiplier>1.0e12) {
2050    freeMultiplier=1.0e12;
2051  } 
2052  freeMultiplier=0.5/freeMultiplier;
2053  double condition = cholesky_->choleskyCondition();
2054  bool caution;
2055  if ((condition<1.0e10&&trueNorm<1.0e12)||numberIterations_<20) {
2056    caution=false;
2057  } else {
2058    caution=true;
2059  } 
2060  // do reduced costs
2061  CoinMemcpyN(dual_,numberRows_,dj_+numberColumns_);
2062  CoinMemcpyN(cost_,numberColumns_,dj_);
2063  matrix_->transposeTimes(-1.0,dual_,dj_);
2064  double extra=eExtra;
2065  double largeGap=1.0e2*solutionNorm_;
2066  if (largeGap<1.0e2) {
2067    largeGap=1.0e2;
2068  } 
2069  double dualFake=0.0;
2070  double dualTolerance =  dblParam_[ClpDualTolerance];
2071  dualTolerance=dualTolerance/scaleFactor_;
2072  if (dualTolerance<1.0e-12) {
2073    dualTolerance=1.0e-12;
2074  } 
2075  double offsetObjective=0.0;
2076  const double killTolerance=primalTolerance();
2077  double qDiagonal;
2078  if (mu_<1.0) {
2079    qDiagonal=1.0e-8;
2080  } else {
2081    qDiagonal=1.0e-8*mu_;
2082  } 
2083  double widenGap=1.0e1;
2084  //largest allowable ratio of lowerSlack/zVec (etc)
2085  double largestRatio;
2086  double epsilonBase;
2087  double diagonalLimit;
2088  if (!caution) {
2089    epsilonBase=eBase;
2090    largestRatio=eRatio;
2091    diagonalLimit=eDiagonal;
2092  } else {
2093    epsilonBase=eBaseCaution;
2094    largestRatio=eRatioCaution;
2095    diagonalLimit=eDiagonalCaution;
2096  } 
2097  double smallGap=1.0e2;
2098  double maximumDJInfeasibility=0.0;
2099  int numberIncreased=0;
2100  int numberDecreased=0;
2101  double largestDiagonal=0.0;
2102  double smallestDiagonal=1.0e50;
2103  double * zVec = zVec_;
2104  double * tVec = tVec_;
2105  double * primal = solution_;
2106  double * dual = dj_;
2107  double * lower = lower_;
2108  double * upper = upper_;
2109  double * lowerSlack = lowerSlack_;
2110  double * upperSlack = upperSlack_;
2111  double * deltaZ = deltaZ_;
2112  double * deltaT = deltaT_;
2113  double * diagonal = diagonal_;
2114  //direction vector in deltaX
2115  double * deltaX = deltaX_;
2116  double * cost = cost_;
2117  for (int iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
2118    if (!flagged(iColumn)) {
2119      double reducedCost=dual[iColumn];
2120      bool thisKilled=false;
2121      double zValue = zVec[iColumn] + actualDualStep_*deltaZ[iColumn];
2122      double tValue = tVec[iColumn] + actualDualStep_*deltaT[iColumn];
2123      zVec[iColumn]=zValue;
2124      tVec[iColumn]=tValue;
2125      double thisWeight=deltaX[iColumn];
2126      double oldPrimal = primal[iColumn];
2127      double newPrimal=primal[iColumn]+actualPrimalStep_*thisWeight;
2128      double primalObjectiveThis=newPrimal*cost[iColumn];
2129      double dualObjectiveThis=0.0;
2130      double t=extra;
2131      double s=extra;
2132      double kill;
2133      if (fabs(newPrimal)>1.0e4) {
2134        kill=killTolerance*1.0e-4*newPrimal;
2135        //kill=killTolerance;
2136        //kill*=1.0e-3;//???????
2137      } else {
2138        kill=killTolerance;
2139        //kill*=1.0e-3;//???????
2140      } 
2141      kill*=1.0e-3;//be conservative
2142      double smallerSlack=COIN_DBL_MAX;
2143      double largerzw;
2144      if (zValue>tValue) {
2145        largerzw=zValue;
2146      } else {
2147        largerzw=tValue;
2148      } 
2149      if (lowerBound(iColumn)) {
2150        double oldSlack = lowerSlack[iColumn];
2151        double newSlack;
2152        if (!fakeLower(iColumn)) {
2153          newSlack=
2154                   lowerSlack[iColumn]+actualPrimalStep_*(oldPrimal-oldSlack
2155                + thisWeight-lower[iColumn]);
2156        } else {
2157          newSlack= lowerSlack[iColumn]+actualPrimalStep_*thisWeight;
2158          if (newSlack<0.0) {
2159            abort();
2160          } 
2161        } 
2162        double epsilon = fabs(newSlack)*epsilonBase;
2163        if (epsilon>1.0e-5) {
2164          //cout<<"bad"<<endl;
2165          epsilon=1.0e-5;
2166        } 
2167        //for changing slack
2168        double zValue2 = zValue;
2169        if (zValue2<epsilonBase) {
2170          zValue2=epsilonBase;
2171        } 
2172        //make sure reasonable
2173        if (zValue<epsilon) {
2174          zValue=epsilon;
2175        } 
2176        //store modified zVec
2177        //no zVec[iColumn]=zValue;
2178        double feasibleSlack=newPrimal-lower[iColumn];
2179        if (feasibleSlack>0.0&&newSlack>0.0) {
2180          double smallGap2=smallGap;
2181          if (fabs(0.1*newPrimal)>smallGap) {
2182            smallGap2=0.1*fabs(newPrimal);
2183          } 
2184          if (!fakeLower(iColumn)) {
2185            double larger;
2186            if (newSlack>feasibleSlack) {
2187              larger=newSlack;
2188            } else {
2189              larger=feasibleSlack;
2190            } 
2191            if (fabs(feasibleSlack-newSlack)<1.0e-6*larger) {
2192              newSlack=feasibleSlack;
2193            } 
2194            //set FAKE here
2195            if (newSlack>zValue2*largestRatio&&newSlack>smallGap2) {
2196              setFakeLower(iColumn);
2197              newSlack=zValue2*largestRatio;
2198              if (newSlack<smallGap2) {
2199                newSlack=smallGap2;
2200              } 
2201              numberDecreased++;
2202            } 
2203          } else {
2204            newSlack=zValue2*largestRatio;
2205            if (newSlack<smallGap2) {
2206              newSlack=smallGap2;
2207            } 
2208            if (newSlack>largeGap) {
2209              //increase up to smaller of z.. and largeGap
2210              newSlack=largeGap;
2211            } 
2212            if (newSlack>widenGap*oldSlack) {
2213              newSlack=widenGap*oldSlack;
2214              numberIncreased++;
2215              //cout<<"wider "<<j<<" "<<newSlack<<" "<<oldSlack<<" "<<feasibleSlack<<endl;
2216            } 
2217            if (newSlack>feasibleSlack) {
2218              newSlack=feasibleSlack;
2219              clearFakeLower(iColumn);
2220            } 
2221          } 
2222        } 
2223        if (zVec[iColumn]>dualTolerance) {
2224          dualObjectiveThis+=lower[iColumn]*zVec[iColumn];
2225        } 
2226        lowerSlack[iColumn]=newSlack;
2227        if (newSlack<smallerSlack) {
2228          smallerSlack=newSlack;
2229        } 
2230        if (!fakeLower(iColumn)) {
2231          double infeasibility = fabs(newPrimal-lowerSlack[iColumn]-lower[iColumn]);
2232          if (infeasibility>maximumBoundInfeasibility) {
2233            maximumBoundInfeasibility=infeasibility;
2234          } 
2235        } 
2236        if (lowerSlack[iColumn]<=kill&&fabs(newPrimal-lower[iColumn])<=kill) {
2237          //may be better to leave at value?
2238          newPrimal=lower[iColumn];
2239          lowerSlack[iColumn]=0.0;
2240          thisKilled=true;
2241          //cout<<j<<" l "<<reducedCost<<" "<<zVec[iColumn]<<endl;
2242        } else {
2243          s+=lowerSlack[iColumn];
2244        } 
2245      } 
2246      if (upperBound(iColumn)) {
2247        //reducedCost-=perturbation;
2248        //primalObjectiveThis-=perturbation*newPrimal;
2249        double oldSlack = upperSlack[iColumn];
2250        double newSlack;
2251        if (!fakeUpper(iColumn)) {
2252          newSlack=
2253                   upperSlack[iColumn]+actualPrimalStep_*(-oldPrimal-oldSlack
2254                - thisWeight+upper[iColumn]);
2255        } else {
2256          newSlack= upperSlack[iColumn]-actualPrimalStep_*thisWeight;
2257        } 
2258        double epsilon = fabs(newSlack)*epsilonBase;
2259        if (epsilon>1.0e-5) {
2260          //cout<<"bad"<<endl;
2261          epsilon=1.0e-5;
2262        } 
2263        //for changing slack
2264        double tValue2 = tValue;
2265        if (tValue2<epsilonBase) {
2266          tValue2=epsilonBase;
2267        } 
2268        //make sure reasonable
2269        if (tValue<epsilon) {
2270          tValue=epsilon;
2271        } 
2272        //store modified tVec
2273        //no tVec[iColumn]=tValue;
2274        double feasibleSlack=upper[iColumn]-newPrimal;
2275        if (feasibleSlack>0.0&&newSlack>0.0) {
2276          double smallGap2=smallGap;
2277          if (fabs(0.1*newPrimal)>smallGap) {
2278            smallGap2=0.1*fabs(newPrimal);
2279          } 
2280          if (!fakeUpper(iColumn)) {
2281            double larger;
2282            if (newSlack>feasibleSlack) {
2283              larger=newSlack;
2284            } else {
2285              larger=feasibleSlack;
2286            } 
2287            if (fabs(feasibleSlack-newSlack)<1.0e-6*larger) {
2288              newSlack=feasibleSlack;
2289            } 
2290            //set FAKE here
2291            if (newSlack>tValue2*largestRatio&&newSlack>smallGap2) {
2292              setFakeUpper(iColumn);
2293              newSlack=tValue2*largestRatio;
2294              if (newSlack<smallGap2) {
2295                newSlack=smallGap2;
2296              } 
2297              numberDecreased++;
2298            } 
2299          } else {
2300            newSlack=tValue2*largestRatio;
2301            if (newSlack<smallGap2) {
2302              newSlack=smallGap2;
2303            } 
2304            if (newSlack>largeGap) {
2305              //increase up to smaller of w.. and largeGap
2306              newSlack=largeGap;
2307            } 
2308            if (newSlack>widenGap*oldSlack) {
2309              numberIncreased++;
2310              //cout<<"wider "<<-j<<" "<<newSlack<<" "<<oldSlack<<" "<<feasibleSlack<<endl;
2311            } 
2312            if (newSlack>feasibleSlack) {
2313              newSlack=feasibleSlack;
2314              clearFakeUpper(iColumn);
2315            } 
2316          } 
2317        } 
2318        if (tVec[iColumn]>dualTolerance) {
2319          dualObjectiveThis-=upper[iColumn]*tVec[iColumn];
2320        } 
2321        upperSlack[iColumn]=newSlack;
2322        if (newSlack<smallerSlack) {
2323          smallerSlack=newSlack;
2324        } 
2325        if (!fakeUpper(iColumn)) {
2326          double infeasibility = fabs(newPrimal+upperSlack[iColumn]-upper[iColumn]);
2327          if (infeasibility>maximumBoundInfeasibility) {
2328            maximumBoundInfeasibility=infeasibility;
2329          } 
2330        } 
2331        if (upperSlack[iColumn]<=kill&&fabs(newPrimal-upper[iColumn])<=kill) {
2332          //may be better to leave at value?
2333          newPrimal=upper[iColumn];
2334          upperSlack[iColumn]=0.0;
2335          thisKilled=true;
2336        } else {
2337          t+=upperSlack[iColumn];
2338        } 
2339      } 
2340      if ((!lowerBound(iColumn))&&
2341               (!upperBound(iColumn))) {
2342        double gap;
2343        if (fabs(newPrimal)>eFree) {
2344          gap=fabs(newPrimal);
2345        } else {
2346          gap=eFree;
2347        } 
2348        zVec[iColumn]=gap*freeMultiplier;
2349        tVec[iColumn]=zVec[iColumn];
2350        //fake to give correct result
2351        s=1.0;
2352        t=1.0;
2353        tValue=0.0;
2354        zValue=2.0*freeMultiplier+qDiagonal;
2355      } 
2356      primal[iColumn]=newPrimal;
2357      if (fabs(newPrimal)>solutionNorm) {
2358        solutionNorm=fabs(newPrimal);
2359      } 
2360      if (!thisKilled) {
2361        primalObjectiveValue+=primalObjectiveThis;
2362        dualObjectiveValue+=dualObjectiveThis;
2363        if (s>largeGap) {
2364          s=largeGap;
2365        } 
2366        if (t>largeGap) {
2367          t=largeGap;
2368        } 
2369        double divisor = s*tValue+t*zValue;
2370        double diagonalValue=(t*s)/divisor;
2371        diagonal[iColumn]=diagonalValue;
2372        //FUDGE
2373        if (diagonalValue>diagonalLimit) {
2374          //cout<<"large diagonal "<<diagonalValue<<endl;
2375          diagonal[iColumn]=diagonalLimit;
2376        } 
2377        if (diagonalValue<1.0e-10) {
2378          //cout<<"small diagonal "<<diagonalValue<<endl;
2379        } 
2380        if (diagonalValue>largestDiagonal) {
2381          largestDiagonal=diagonalValue;
2382        } 
2383        if (diagonalValue<smallestDiagonal) {
2384          smallestDiagonal=diagonalValue;
2385        } 
2386        double dualInfeasibility=reducedCost-zVec[iColumn]+tVec[iColumn];
2387        if (fabs(dualInfeasibility)>dualTolerance) {
2388          if ((!lowerBound(iColumn))&&
2389                 (!upperBound(iColumn))) {
2390            dualFake+=newPrimal*dualInfeasibility;
2391          } else {
2392            dualFake+=newPrimal*dualInfeasibility;
2393          } 
2394        } 
2395        dualInfeasibility=fabs(dualInfeasibility);
2396        if (dualInfeasibility>maximumDualError) {
2397          maximumDualError=dualInfeasibility;
2398        } 
2399        deltaZ[iColumn]=0.0;
2400      } else {
2401        numberKilled++;
2402        diagonal[iColumn]=0.0;
2403        zVec[iColumn]=0.0;
2404        tVec[iColumn]=0.0;
2405        setFlagged(iColumn);
2406        setFixedOrFree(iColumn);
2407        deltaZ[iColumn]=newPrimal;
2408        offsetObjective+=newPrimal*cost[iColumn];
2409      } 
2410    } else {
2411      deltaZ[iColumn]=primal[iColumn];
2412      diagonal[iColumn]=0.0;
2413      offsetObjective+=primal[iColumn]*cost[iColumn];
2414      if (upper[iColumn]-lower[iColumn]>1.0e-5) {
2415        if (primal[iColumn]<lower[iColumn]+1.0e-8&&dual[iColumn]<-1.0e-8) {
2416          if (-dual[iColumn]>maximumDJInfeasibility) {
2417            maximumDJInfeasibility=-dual[iColumn];
2418          } 
2419        } 
2420        if (primal[iColumn]>upper[iColumn]-1.0e-8&&dual[iColumn]>1.0e-8) {
2421          if (dual[iColumn]>maximumDJInfeasibility) {
2422            maximumDJInfeasibility=dual[iColumn];
2423          } 
2424        } 
2425      } 
2426    } 
2427  } 
2428  handler_->message(CLP_BARRIER_DIAGONAL,messages_)
2429    <<largestDiagonal<<smallestDiagonal
2430    <<CoinMessageEol;
2431#if 0
2432  // If diagonal wild - kill some
2433  if (largestDiagonal>1.0e17*smallestDiagonal) {
2434    double killValue =largestDiagonal*1.0e-17;
2435    for (int iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
2436      if (fabs(diagonal_[iColumn])<killValue)
2437        diagonal_[iolumn]=0.0;
2438    }
2439  }
2440#endif
2441  // update rhs region
2442  multiplyAdd(deltaZ+numberColumns_,numberRows_,-1.0,rhsFixRegion_,1.0);
2443  matrix_->times(1.0,deltaZ,rhsFixRegion_);
2444  primalObjectiveValue=primalObjectiveValue+offsetObjective;
2445  dualObjectiveValue+=offsetObjective+dualFake;
2446  if (numberIncreased||numberDecreased) {
2447    handler_->message(CLP_BARRIER_SLACKS,messages_)
2448      <<numberIncreased<<numberDecreased
2449      <<CoinMessageEol;
2450  } 
2451  if (maximumDJInfeasibility) {
2452    handler_->message(CLP_BARRIER_DUALINF,messages_)
2453      <<maximumDJInfeasibility
2454      <<CoinMessageEol;
2455  } 
2456  // Need to rethink (but it is only for printing)
2457  sumPrimalInfeasibilities_=maximumRhsInfeasibility;
2458  sumDualInfeasibilities_=maximumDualError;
2459  maximumBoundInfeasibility_ = maximumBoundInfeasibility;
2460  //compute error and fixed RHS
2461  multiplyAdd(solution_+numberColumns_,numberRows_,-1.0,errorRegion_,0.0);
2462  matrix_->times(1.0,solution_,errorRegion_);
2463  maximumDualError_=maximumDualError;
2464  maximumBoundInfeasibility_=maximumBoundInfeasibility;
2465  solutionNorm_=solutionNorm;
2466  //finish off objective computation
2467  primalObjective_=primalObjectiveValue*scaleFactor_;
2468  double dualValue2=innerProduct(dual_,numberRows_,
2469                rhsFixRegion_);
2470  dualObjectiveValue-=dualValue2;
2471  dualObjective_=dualObjectiveValue*scaleFactor_;
2472  if (numberKilled) {
2473      handler_->message(CLP_BARRIER_KILLED,messages_)
2474        <<numberKilled
2475        <<CoinMessageEol;
2476  } 
2477  double maximumRHSError1=0.0;
2478  double maximumRHSError2=0.0;
2479  double primalOffset=0.0;
2480  char * dropped = cholesky_->rowsDropped();
2481  int iRow;
2482  for (iRow=0;iRow<numberRows_;iRow++) {
2483    double value=errorRegion_[iRow];
2484    if (!dropped[iRow]) {
2485      if (fabs(value)>maximumRHSError1) {
2486        maximumRHSError1=fabs(value);
2487      } 
2488    } else {
2489      if (fabs(value)>maximumRHSError2) {
2490        maximumRHSError2=fabs(value);
2491      } 
2492      primalOffset+=value*dual_[iRow];
2493    } 
2494  } 
2495  primalObjective_-=primalOffset*scaleFactor_;
2496  if (maximumRHSError1>maximumRHSError2) {
2497    maximumRHSError_=maximumRHSError1;
2498  } else {
2499    maximumRHSError_=maximumRHSError1; //note change
2500    if (maximumRHSError2>primalTolerance()) {
2501      handler_->message(CLP_BARRIER_ABS_DROPPED,messages_)
2502        <<maximumRHSError2
2503        <<CoinMessageEol;
2504    } 
2505  } 
2506  objectiveNorm_=maximumAbsElement(dual_,numberRows_);
2507  if (objectiveNorm_<1.0e-12) {
2508    objectiveNorm_=1.0e-12;
2509  } 
2510  if (objectiveNorm_<baseObjectiveNorm_) {
2511    //std::cout<<" base "<<baseObjectiveNorm_<<" "<<objectiveNorm_<<std::endl;
2512    if (objectiveNorm_<baseObjectiveNorm_*1.0e-4) {
2513      objectiveNorm_=baseObjectiveNorm_*1.0e-4;
2514    } 
2515  } 
2516  bool primalFeasible=true;
2517  if (maximumRHSError_>primalTolerance()||
2518      maximumDualError_>dualTolerance/scaleFactor_) {
2519    handler_->message(CLP_BARRIER_ABS_ERROR,messages_)
2520      <<maximumRHSError_<<maximumDualError_
2521      <<CoinMessageEol;
2522  } 
2523  if (rhsNorm_>solutionNorm_) {
2524    solutionNorm_=rhsNorm_;
2525  } 
2526  double scaledRHSError=maximumRHSError_/solutionNorm_;
2527  bool dualFeasible=true;
2528  if (maximumBoundInfeasibility_>primalTolerance()||
2529      scaledRHSError>primalTolerance())
2530    primalFeasible=false;
2531  if (maximumDualError_>objectiveNorm_*dualTolerance) 
2532    dualFeasible=false;
2533  if (!primalFeasible||!dualFeasible) {
2534    handler_->message(CLP_BARRIER_FEASIBLE,messages_)
2535      <<maximumBoundInfeasibility_<<scaledRHSError
2536      <<maximumDualError_/objectiveNorm_
2537      <<CoinMessageEol;
2538  } 
2539  if (!gonePrimalFeasible_) {
2540    gonePrimalFeasible_=primalFeasible;
2541  } else if (!primalFeasible) {
2542    gonePrimalFeasible_=primalFeasible;
2543    if (!numberKilled) {
2544      handler_->message(CLP_BARRIER_GONE_INFEASIBLE,messages_)
2545      <<CoinMessageEol;
2546    } 
2547  } 
2548  if (!goneDualFeasible_) {
2549    goneDualFeasible_=dualFeasible;
2550  } else if (!dualFeasible) {
2551    handler_->message(CLP_BARRIER_GONE_INFEASIBLE,messages_)
2552      <<CoinMessageEol;
2553    goneDualFeasible_=dualFeasible;
2554  } 
2555  //objectiveValue();
2556  if (solutionNorm_>1.0e40) {
2557    std::cout <<"primal off to infinity"<<std::endl;
2558    abort();
2559  } 
2560  if (objectiveNorm_>1.0e40) {
2561    std::cout <<"dual off to infinity"<<std::endl;
2562    abort();
2563  } 
2564  handler_->message(CLP_BARRIER_STEP,messages_)
2565    <<actualPrimalStep_
2566    <<actualDualStep_
2567    <<mu_
2568    <<CoinMessageEol;
2569  numberIterations_++;
2570  return numberKilled;
2571}
2572//  Save info on products of affine deltaT*deltaW and deltaS*deltaZ
2573double 
2574ClpPredictorCorrector::affineProduct()
2575{
2576  double * deltaZ = deltaZ_;
2577  double * deltaT = deltaT_;
2578  double * lower = lower_;
2579  double * upper = upper_;
2580  double * lowerSlack = lowerSlack_;
2581  double * upperSlack = upperSlack_;
2582  double * primal = solution_;
2583  //direction vector in deltaX
2584  double * deltaX = deltaX_;
2585  double product = 0.0;
2586  //IF zVec starts as 0 then deltaZ always zero
2587  //(remember if free then zVec not 0)
2588  //I think free can be done with careful use of boundSlacks to zero
2589  //out all we want
2590  for (int iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
2591    double w3=deltaZ[iColumn]*deltaX[iColumn];
2592    double w4=-deltaT[iColumn]*deltaX[iColumn];
2593    if (lowerBound(iColumn)) {
2594      if (!fakeLower(iColumn)) {
2595        w3+=deltaZ[iColumn]*(primal[iColumn]-lowerSlack[iColumn]-lower[iColumn]);
2596      } 
2597      product+=w3;
2598    } 
2599    if (upperBound(iColumn)) {
2600      if (!fakeUpper(iColumn)) {
2601        w4+=deltaT[iColumn]*(-primal[iColumn]-upperSlack[iColumn]+upper[iColumn]);
2602      } 
2603      product+=w4;
2604    } 
2605  } 
2606  return product;
2607}
Note: See TracBrowser for help on using the repository browser.