source: trunk/ClpPredictorCorrector.cpp @ 365

Last change on this file since 365 was 365, checked in by jpfasano, 16 years ago

Modified to compile with MS VC++ Ver 6

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