CLP: ClpPredictorCorrector.cpp Source File

00001 // Copyright (C) 2003, International Business Machines
00002 // Corporation and others.  All Rights Reserved.
00003 
00004 /* 
00005    Implements crude primal dual predictor corrector algorithm
00006 
00007  */
00008 
00009 
00010 #include "CoinPragma.hpp"
00011 
00012 #include <math.h>
00013 
00014 #include "CoinHelperFunctions.hpp"
00015 #include "ClpPredictorCorrector.hpp"
00016 #include "CoinPackedMatrix.hpp"
00017 #include "ClpMessage.hpp"
00018 #include "ClpCholeskyDense.hpp"
00019 #include "ClpHelperFunctions.hpp"
00020 #include <cfloat>
00021 #include <cassert>
00022 #include <string>
00023 #include <cstdio>
00024 #include <iostream>
00025 
00026 static double eScale=1.0e57;
00027 static double eBaseCaution=1.0e-12;
00028 static double eBase=1.0e-12;
00029 static double eRatio=1.0e40;
00030 static double eRatioCaution=1.0e25;
00031 static double eDiagonal=1.0e25;
00032 static double eDiagonalCaution=1.0e18;
00033 static double eExtra=1.0e-12;
00034 static double eFree =1.0e3;
00035 
00036 // main function
00037 
00038 int ClpPredictorCorrector::solve ( )
00039 {
00040   problemStatus_=-1;
00041   algorithm_=1;
00042   //create all regions
00043   if (!createWorkingData()) {
00044     problemStatus_=4;
00045     return 2;
00046   }
00047   //initializeFeasible(); - this just set fixed flag
00048   smallestInfeasibility_=COIN_DBL_MAX;
00049   int i;
00050   for (i=0;i<LENGTH_HISTORY;i++) 
00051     historyInfeasibility_[i]=COIN_DBL_MAX;
00052 
00053   //bool firstTime=true;
00054   //firstFactorization(true);
00055   cholesky_->order(this);
00056   mu_=1.0e10;
00057   //set iterations
00058   numberIterations_=-1;
00059   //initialize solution here
00060   createSolution();
00061   //firstFactorization(false);
00062   CoinZeroN(dual_,numberRows_);
00063   multiplyAdd(solution_+numberColumns_,numberRows_,-1.0,errorRegion_,0.0);
00064   matrix_->times(1.0,solution_,errorRegion_);
00065   maximumRHSError_=maximumAbsElement(errorRegion_,numberRows_);
00066   maximumBoundInfeasibility_=maximumRHSError_;
00067   //double maximumDualError_=COIN_DBL_MAX;
00068   //initialize
00069   actualDualStep_=0.0;
00070   actualPrimalStep_=0.0;
00071   gonePrimalFeasible_=false;
00072   goneDualFeasible_=false;
00073   //bool hadGoodSolution=false;
00074   diagonalScaleFactor_=1.0;
00075   diagonalNorm_=solutionNorm_;
00076   mu_=solutionNorm_;
00077   int numberFixed=updateSolution();
00078   int numberFixedTotal=numberFixed;
00079   //int numberRows_DroppedBefore=0;
00080   //double extra=eExtra;
00081   //double maximumPerturbation=COIN_DBL_MAX;
00082   //constants for infeas interior point
00083   const double beta2 = 0.99995;
00084   const double tau   = 0.00002;
00085   double lastComplementarityGap=COIN_DBL_MAX;
00086   int lastGoodIteration=0;
00087   double bestObjectiveGap=COIN_DBL_MAX;
00088   int saveIteration=-1;
00089   bool sloppyOptimal=false;
00090   double * savePi=NULL;
00091   double * savePrimal=NULL;
00092   int numberTotal = numberRows_+numberColumns_;
00093   while (problemStatus_<0) {
00094     complementarityGap_=complementarityGap(numberComplementarityPairs_,0);
00095       handler_->message(CLP_BARRIER_ITERATION,messages_)
00096     <<numberIterations_
00097     <<primalObjective_*optimizationDirection_- dblParam_[ClpObjOffset]
00098     << dualObjective_*optimizationDirection_- dblParam_[ClpObjOffset]
00099     <<complementarityGap_
00100     <<numberFixedTotal
00101     <<cholesky_->rank()
00102     <<CoinMessageEol;
00103     double goodGapChange;
00104     if (!sloppyOptimal) {
00105       goodGapChange=0.93;
00106     } else {
00107       goodGapChange=0.7;
00108     } 
00109     double gapO;
00110     double lastGood=bestObjectiveGap;
00111     if (gonePrimalFeasible_&&goneDualFeasible_) {
00112       double largestObjective;
00113       if (fabs(primalObjective_)>fabs(dualObjective_)) {
00114         largestObjective = fabs(primalObjective_);
00115       } else {
00116         largestObjective = fabs(dualObjective_);
00117       } 
00118       if (largestObjective<1.0) {
00119         largestObjective=1.0;
00120       } 
00121       gapO=fabs(primalObjective_-dualObjective_)/largestObjective;
00122       handler_->message(CLP_BARRIER_OBJECTIVE_GAP,messages_)
00123         <<gapO
00124         <<CoinMessageEol;
00125       //start saving best
00126       if (gapO<bestObjectiveGap) {
00127         saveIteration=numberIterations_;
00128         bestObjectiveGap=gapO;
00129         if (!savePi) {
00130           savePi=new double[numberRows_];
00131           savePrimal = new double [numberTotal];
00132         } 
00133         CoinMemcpyN(dual_,numberRows_,savePi);
00134         CoinMemcpyN(solution_,numberTotal,savePrimal);
00135       } else if(gapO>2.0*bestObjectiveGap) {
00136         //maybe be more sophisticated e.g. re-initialize having
00137         //fixed variables and dropped rows
00138         //std::cout <<" gap increasing "<<std::endl;
00139       } 
00140       //std::cout <<"could stop"<<std::endl;
00141       //gapO=0.0;
00142       if (fabs(primalObjective_-dualObjective_)<dualTolerance()) {
00143         gapO=0.0;
00144       } 
00145     } else {
00146       gapO=COIN_DBL_MAX;
00147       if (saveIteration>=0) {
00148         handler_->message(CLP_BARRIER_GONE_INFEASIBLE,messages_)
00149           <<CoinMessageEol;
00150       } 
00151     } 
00152     if (gapO<1.0e-7&&!sloppyOptimal) {
00153       sloppyOptimal=true;
00154       handler_->message(CLP_BARRIER_CLOSE_TO_OPTIMAL,messages_)
00155         <<numberIterations_<<complementarityGap_
00156         <<CoinMessageEol;
00157     } 
00158     if (complementarityGap_>=1.05*lastComplementarityGap) {
00159       handler_->message(CLP_BARRIER_COMPLEMENTARITY,messages_)
00160         <<complementarityGap_<<"increasing"
00161         <<CoinMessageEol;
00162       if (saveIteration>=0&&sloppyOptimal) {
00163         handler_->message(CLP_BARRIER_EXIT2,messages_)
00164           <<saveIteration
00165           <<CoinMessageEol;
00166         break;
00167       } else {
00168         //lastComplementarityGap=complementarityGap_;
00169       } 
00170     } else if (complementarityGap_<goodGapChange*lastComplementarityGap) {
00171       lastGoodIteration=numberIterations_;
00172       lastComplementarityGap=complementarityGap_;
00173     } else if (numberIterations_-lastGoodIteration>=5&&complementarityGap_<1.0e-3) {
00174       handler_->message(CLP_BARRIER_COMPLEMENTARITY,messages_)
00175         <<complementarityGap_<<"not decreasing"
00176         <<CoinMessageEol;
00177       if (gapO>0.75*lastGood) {
00178         break;
00179       } 
00180     } 
00181     if (numberIterations_>maximumBarrierIterations_) {
00182       handler_->message(CLP_BARRIER_STOPPING,messages_)
00183         <<CoinMessageEol;
00184       break;
00185     } 
00186     if (gapO<targetGap_) {
00187       problemStatus_=0;
00188       handler_->message(CLP_BARRIER_EXIT,messages_)
00189         <<" "
00190         <<CoinMessageEol;
00191         break;//finished
00192     } 
00193     if (complementarityGap_<1.0e-18) {
00194       problemStatus_=0;
00195       handler_->message(CLP_BARRIER_EXIT,messages_)
00196         <<"- small complementarity gap"
00197         <<CoinMessageEol;
00198         break;//finished
00199     } 
00200     if (complementarityGap_<1.0e-10&&gapO<1.0e-10) {
00201       problemStatus_=0;
00202       handler_->message(CLP_BARRIER_EXIT,messages_)
00203         <<"- objective gap and complementarity gap both small"
00204         <<CoinMessageEol;
00205         break;//finished
00206     } 
00207     if (gapO<1.0e-9) {
00208       double value=gapO*complementarityGap_;
00209       value*=actualPrimalStep_;
00210       value*=actualDualStep_;
00211       //std::cout<<value<<std::endl;
00212       if (value<1.0e-17&&numberIterations_>lastGoodIteration) {
00213         problemStatus_=0;
00214         handler_->message(CLP_BARRIER_EXIT,messages_)
00215           <<"- objective gap and complementarity gap both smallish and small steps"
00216           <<CoinMessageEol;
00217         break;//finished
00218       } 
00219     } 
00220     bool useAffine=false;
00221     bool goodMove=false;
00222     bool doCorrector=true;
00223     //bool retry=false;
00224     worstDirectionAccuracy_=0.0;
00225     while (!goodMove) {
00226       goodMove=true;
00227       int newDropped=0;
00228       //Predictor step
00229       //Are we going to use the affine direction?
00230       if (!useAffine) {
00231         //no - normal
00232         //prepare for cholesky.  Set up scaled diagonal in weights
00233         //  ** for efficiency may be better if scale factor known before
00234         double norm2=0.0;
00235         double maximumValue;
00236         getNorms(diagonal_,numberColumns_,maximumValue,norm2);
00237         diagonalNorm_ = sqrt(norm2/numberComplementarityPairs_);
00238         diagonalScaleFactor_=1.0;
00239         double maximumAllowable=eScale;
00240         //scale so largest is less than allowable ? could do better
00241         double factor=0.5;
00242         while (maximumValue>maximumAllowable) {
00243           diagonalScaleFactor_*=factor;
00244           maximumValue*=factor;
00245         } /* endwhile */
00246         if (diagonalScaleFactor_!=1.0) {
00247           handler_->message(CLP_BARRIER_SCALING,messages_)
00248             <<"diagonal"<<diagonalScaleFactor_
00249             <<CoinMessageEol;
00250           diagonalNorm_*=diagonalScaleFactor_;
00251         } 
00252         multiplyAdd(NULL,numberTotal,0.0,diagonal_,
00253                       diagonalScaleFactor_);
00254         int * rowsDroppedThisTime = new int [numberRows_];
00255         newDropped=cholesky_->factorize(diagonal_,rowsDroppedThisTime);
00256         if (newDropped) {
00257           //int newDropped2=cholesky_->factorize(diagonal_,rowsDroppedThisTime);
00258           //assert(!newDropped2);
00259           if (newDropped<0) {
00260             //replace dropped
00261             newDropped=-newDropped;
00262             //off 1 to allow for reset all
00263             newDropped--;
00264             //set all bits false
00265             cholesky_->resetRowsDropped();
00266           } 
00267         } 
00268         delete [] rowsDroppedThisTime;
00269         if (cholesky_->status()) {
00270           std::cout << "bad cholesky?" <<std::endl;
00271           abort();
00272         } 
00273       } 
00274       //set up for affine direction
00275       setupForSolve(0);
00276       double directionAccuracy=findDirectionVector(0);
00277       if (directionAccuracy>worstDirectionAccuracy_) {
00278         worstDirectionAccuracy_=directionAccuracy;
00279       } 
00280       int phase=0; // predictor, corrector , primal dual
00281       // 0 - normal
00282       // 1 - second time around i.e. no need for first part
00283       // 2 - affine step only
00284       // 9 - to exit from while because of error (to go round again)
00285       // (9 also used to signal end of while (but then goodMove is true))
00286       int recoveryMode=0;
00287       if (!goodMove) {
00288         recoveryMode=9;
00289       } 
00290       if (goodMove&&useAffine) {
00291         recoveryMode=2;
00292         phase=0;
00293       } 
00294       while (recoveryMode!=9) {
00295         goodMove=true;
00296         if (!recoveryMode) {
00297           findStepLength(phase);
00298           int nextNumber; //number of complementarity pairs
00299           double nextGap=complementarityGap(nextNumber,1);
00300           //if (complementarityGap_>1.0e-1&&worstDirectionAccuracy_<1.0e-5) {
00301           //wasif (complementarityGap_>1.0e-2*numberComplementarityPairs_) {
00302           if (complementarityGap_>1.0e-4*numberComplementarityPairs_) {
00303             //std::cout <<"predicted duality gap "<<nextGap<<std::endl;
00304             //double part1=nextGap/numberComplementarityPairs_;
00305             double part1=nextGap/numberComplementarityPairs_;
00306             double part2=nextGap/complementarityGap_;
00307             mu_=part1*part2*part2;
00308           } else {
00309             double phi;
00310             if (numberComplementarityPairs_<=500) {
00311               phi=pow(numberComplementarityPairs_,2);
00312             } else {
00313               phi=pow(numberComplementarityPairs_,1.5);
00314               if (phi<500.0*500.0) {
00315                 phi=500.0*500.0;
00316               } 
00317             }
00318             mu_=complementarityGap_/phi;
00319             //? could use gapO?
00320             //if small then should be stopping
00321             //if (mu_<1.0e-4/(numberComplementarityPairs_*qqqq)) {
00322             //mu_=1.0e-4/(numberComplementarityPairs_*qqqq);
00323             //? better to skip corrector?
00324             //} 
00325           } 
00326           //save information
00327           double product=affineProduct();
00328           //can we do corrector step?
00329           double xx= complementarityGap_*(beta2-tau) +product;
00330           //std::cout<<"gap part "<<
00331           //complementarityGap_*(beta2-tau)<<" after adding product = "<<xx;
00332           if (xx>0.0) {
00333             double saveMu = mu_;
00334             double mu2=numberComplementarityPairs_;
00335             mu2=xx/mu2;
00336             if (mu2>mu_) {
00337               //std::cout<<" could increase to "<<mu2<<std::endl;
00338               //was mu2=mu2*0.25;
00339               mu2=mu2*0.99;
00340               if (mu2<mu_) {
00341                 mu_=mu2;
00342                 //std::cout<<" changing mu to "<<mu_<<std::endl;
00343               } else {
00344                 //std::cout<<std::endl;
00345               } 
00346             } else {
00347               //std::cout<<" should decrease to "<<mu2<<std::endl;
00348               mu_=0.5*mu2;
00349               std::cout<<" changing mu to "<<mu_<<std::endl;
00350             } 
00351             handler_->message(CLP_BARRIER_MU,messages_)
00352               <<saveMu<<mu_
00353               <<CoinMessageEol;
00354           } else {
00355             //std::cout<<" bad by any standards"<<std::endl;
00356           } 
00357           if (complementarityGap_*(beta2-tau)+product-mu_*numberComplementarityPairs_<0.0) {
00358             doCorrector=false;
00359             double floatNumber = numberComplementarityPairs_;
00360             mu_=complementarityGap_/(floatNumber*floatNumber);
00361             //? if small we should be stopping
00362             //if (mu_<1.0e-4/(numberComplementarityPairs_*numberComplementarityPairs_)) {
00363             //mu_=1.0e-4/(totalVariables*totalVariables);
00364             //} 
00365             handler_->message(CLP_BARRIER_INFO,messages_)
00366               <<"no corrector step"
00367               <<CoinMessageEol;
00368             phase=2;
00369           } else {
00370             phase=1;
00371           } 
00372         } 
00373         //set up for next step
00374         setupForSolve(phase);
00375         double directionAccuracy2=findDirectionVector(phase);
00376         if (directionAccuracy2>worstDirectionAccuracy_) {
00377           worstDirectionAccuracy_=directionAccuracy2;
00378         } 
00379         double testValue=1.0e2*directionAccuracy;
00380         if (1.0e2*projectionTolerance_>testValue) {
00381           testValue=1.0e2*projectionTolerance_;
00382         } 
00383         if (primalTolerance()>testValue) {
00384           testValue=primalTolerance();
00385         } 
00386         if (maximumRHSError_>testValue) {
00387           testValue=maximumRHSError_;
00388         } 
00389         if (!recoveryMode) {
00390           if (directionAccuracy2>testValue&&numberIterations_>=-77) {
00391             goodMove=false;
00392             useAffine=true;//if bad accuracy
00393             doCorrector=false;
00394             recoveryMode=9;
00395           } 
00396         } 
00397         if (goodMove) {
00398           findStepLength(phase);
00399           if (numberIterations_>=-77) {
00400             goodMove=checkGoodMove(doCorrector);
00401           } else {
00402             goodMove=true;
00403           } 
00404           if (!goodMove) {
00405             if (doCorrector) {
00406               doCorrector=false;
00407               double floatNumber = numberComplementarityPairs_;
00408               mu_=complementarityGap_/(floatNumber*floatNumber);
00409               handler_->message(CLP_BARRIER_INFO,messages_)
00410                 <<" no corrector step - original move would be bad"
00411                 <<CoinMessageEol;
00412               phase=2;
00413               recoveryMode=1;
00414             } else {
00415               // if any killed then do zero step and hope for best
00416               abort();
00417             } 
00418           } 
00419         } 
00420         //force leave
00421         if (goodMove) {
00422           recoveryMode=9;
00423         } 
00424       } /* endwhile */
00425     } /* endwhile */
00426     numberFixed=updateSolution();
00427     numberFixedTotal+=numberFixed;
00428   } /* endwhile */
00429   if (savePi) {
00430     //std::cout<<"Restoring from iteration "<<saveIteration<<std::endl;
00431     CoinMemcpyN(savePi,numberRows_,dual_);
00432     CoinMemcpyN(savePrimal,numberTotal,solution_);
00433     delete [] savePi;
00434     delete [] savePrimal;
00435   } 
00436   //recompute slacks
00437   // Split out solution
00438   CoinZeroN(rowActivity_,numberRows_);
00439   CoinMemcpyN(solution_,numberColumns_,columnActivity_);
00440   matrix_->times(1.0,columnActivity_,rowActivity_);
00441   //unscale objective
00442   multiplyAdd(NULL,numberTotal,0.0,cost_,scaleFactor_);
00443   multiplyAdd(NULL,numberRows_,0,dual_,scaleFactor_);
00444   CoinMemcpyN(cost_,numberColumns_,reducedCost_);
00445   matrix_->transposeTimes(-1.0,dual_,reducedCost_);
00446   CoinMemcpyN(reducedCost_,numberColumns_,dj_);
00447   checkSolution();
00448   handler_->message(CLP_BARRIER_END,messages_)
00449     <<sumPrimalInfeasibilities_
00450     <<sumDualInfeasibilities_
00451     <<complementarityGap_
00452     <<objectiveValue()
00453     <<CoinMessageEol;
00454   //std::cout<<"Absolute primal infeasibility at end "<<sumPrimalInfeasibilities_<<std::endl;
00455   //std::cout<<"Absolute dual infeasibility at end "<<sumDualInfeasibilities_<<std::endl;
00456   //std::cout<<"Absolute complementarity at end "<<complementarityGap_<<std::endl;
00457   //std::cout<<"Primal objective "<<objectiveValue()<<std::endl;
00458   //std::cout<<"maximum complementarity "<<worstComplementarity_<<std::endl;
00459   //delete all temporary regions
00460   deleteWorkingData();
00461   return problemStatus_;
00462 }
00463 // findStepLength.
00464 //phase  - 0 predictor
00465 //         1 corrector
00466 //         2 primal dual
00467 double ClpPredictorCorrector::findStepLength(const int phase)
00468 {
00469   double directionNorm=0.0;
00470   double maximumPrimalStep=COIN_DBL_MAX;
00471   double maximumDualStep=COIN_DBL_MAX;
00472   int numberTotal = numberRows_+numberColumns_;
00473   double tolerance = 1.0e-12;
00474   int chosenPrimalSequence=-1;
00475   int chosenDualSequence=-1;
00476   double extra=eExtra;
00477   double * zVec = zVec_;
00478   double * wVec = wVec_;
00479   double * primal = solution_;
00480   double * lower = lower_;
00481   double * upper = upper_;
00482   double * lowerSlack = lowerSlack_;
00483   double * upperSlack = upperSlack_;
00484   double * work1 = deltaZ_;
00485   double * work2 = deltaW_;
00486   double * work3 = deltaS_;
00487   double * work4 = deltaT_;
00488   //direction vector in weights
00489   double * weights = weights_;
00490   if (!phase) {
00491     //Now get affine deltas for Z(duals on LBds) and W (duals on UBds)
00492     for (int iColumn=0;iColumn<numberTotal;iColumn++) {
00493       if (!flagged(iColumn)) {
00494         double z1=-zVec[iColumn];
00495         double w1=-wVec[iColumn];
00496         double work3Value=0.0;
00497         double work4Value=0.0;
00498         double directionElement=weights[iColumn];
00499         double value=primal[iColumn];
00500         if (directionNorm<fabs(directionElement)) {
00501           directionNorm=fabs(directionElement);
00502         } 
00503         //below does not feel right - can it be simplified because
00504         //of zero values for zVec and wVec
00505         if (lowerBound(iColumn)||
00506                         upperBound(iColumn)) {
00507           if (lowerBound(iColumn)) {
00508             double gap=lowerSlack[iColumn]+extra;
00509             double delta;
00510             if (!fakeLower(iColumn)) {
00511               delta = lower[iColumn]+lowerSlack[iColumn]
00512              -  value - directionElement;
00513             } else {
00514               delta = - directionElement;
00515             } 
00516             z1+=(zVec[iColumn]*delta)/gap;
00517             work3Value=-delta;
00518             if (lowerSlack[iColumn]<maximumPrimalStep*delta) {
00519               maximumPrimalStep=lowerSlack[iColumn]/delta;
00520               chosenPrimalSequence=iColumn;
00521             } 
00522             if (zVec[iColumn]>tolerance) {
00523               if (zVec[iColumn]<-z1*maximumDualStep) {
00524                 maximumDualStep=-zVec[iColumn]/z1;
00525                 chosenDualSequence=iColumn;
00526               } 
00527             } 
00528           } 
00529           if (upperBound(iColumn)) {
00530             double gap=upperSlack[iColumn]+extra;
00531             double delta;
00532             if (!fakeUpper(iColumn)) {
00533               delta = -upper[iColumn]+upperSlack[iColumn]
00534                         + value + directionElement;
00535             } else {
00536               delta = directionElement;
00537             } 
00538             w1+=(wVec[iColumn]*delta)/gap;
00539             work4Value=-delta;
00540             if (upperSlack[iColumn]<maximumPrimalStep*delta) {
00541               maximumPrimalStep=upperSlack[iColumn]/delta;
00542               chosenPrimalSequence=iColumn;
00543             } 
00544             if (wVec[iColumn]>tolerance) {
00545               if (wVec[iColumn]<-w1*maximumDualStep) {
00546                 maximumDualStep=-wVec[iColumn]/w1;
00547                 chosenDualSequence=iColumn;
00548               } 
00549             } 
00550           } 
00551         } else {
00552           //free
00553           double gap=fabs(value);
00554           double multiplier=1.0/gap;
00555           if (gap<1.0) {
00556             multiplier=1,0;
00557           } 
00558           z1-=multiplier*directionElement*zVec[iColumn];
00559           w1+=multiplier*directionElement*wVec[iColumn];
00560         } 
00561         work1[iColumn]=z1;
00562         work2[iColumn]=w1;
00563         work3[iColumn]=work3Value;
00564         work4[iColumn]=work4Value;
00565       } else {
00566         work1[iColumn]=0.0;
00567         work2[iColumn]=0.0;
00568         work3[iColumn]=0.0;
00569         work4[iColumn]=0.0;
00570       } 
00571     } 
00572   } else if (phase==1) {
00573     //corrector step
00574     for (int iColumn=0;iColumn<numberTotal;iColumn++) {
00575       if (!flagged(iColumn)) {
00576         double z1=-zVec[iColumn];
00577         double w1=-wVec[iColumn];
00578         double work3Value=0.0;
00579         double work4Value=0.0;
00580         double directionElement=weights[iColumn];
00581         double value=primal[iColumn];
00582         if (directionNorm<fabs(directionElement)) {
00583           directionNorm=fabs(directionElement);
00584         } 
00585         //below does not feel right - can it be simplified because
00586         //of zero values for zVec and wVec
00587         if (lowerBound(iColumn)||
00588                         upperBound(iColumn)) {
00589           if (lowerBound(iColumn)) {
00590             double gap=lowerSlack[iColumn]+extra;
00591             double delta;
00592             if (!fakeLower(iColumn)) {
00593               delta = lower[iColumn]+lowerSlack[iColumn]
00594              -  value - directionElement;
00595             } else {
00596               delta = - directionElement;
00597             } 
00598             z1+=(mu_-work3[iColumn]+zVec[iColumn]*delta)/gap;
00599             work3Value=-delta;
00600             if (lowerSlack[iColumn]<maximumPrimalStep*delta) {
00601               maximumPrimalStep=lowerSlack[iColumn]/delta;
00602               chosenPrimalSequence=iColumn;
00603             } 
00604             if (zVec[iColumn]>tolerance) {
00605               if (zVec[iColumn]<-z1*maximumDualStep) {
00606                 maximumDualStep=-zVec[iColumn]/z1;
00607                 chosenDualSequence=iColumn;
00608               } 
00609             } 
00610           } 
00611           if (upperBound(iColumn)) {
00612             double gap=upperSlack[iColumn]+extra;
00613             double delta;
00614             if (!fakeUpper(iColumn)) {
00615               delta = -upper[iColumn]+upperSlack[iColumn]
00616                         + value + directionElement;
00617             } else {
00618               delta = directionElement;
00619             } 
00620             //double delta = -upper[iColumn]+upperSlack[iColumn]-
00621               //+ value + directionElement;
00622             w1+=(mu_-work4[iColumn]+wVec[iColumn]*delta)/gap;
00623             work4Value=-delta;
00624             if (upperSlack[iColumn]<maximumPrimalStep*delta) {
00625               maximumPrimalStep=upperSlack[iColumn]/delta;
00626               chosenPrimalSequence=iColumn;
00627             } 
00628             if (wVec[iColumn]>tolerance) {
00629               if (wVec[iColumn]<-w1*maximumDualStep) {
00630                 maximumDualStep=-wVec[iColumn]/w1;
00631                 chosenDualSequence=iColumn;
00632               } 
00633             } 
00634           } 
00635         } else {
00636           //free
00637           double gap=fabs(value);
00638           double multiplier=1.0/gap;
00639           if (gap<1.0) {
00640             multiplier=1,0;
00641           } 
00642           z1+=multiplier*(mu_-work3[iColumn]-directionElement*zVec[iColumn]);
00643           w1+=multiplier*(mu_-work4[iColumn]+directionElement*wVec[iColumn]);
00644         } 
00645         work1[iColumn]=z1;
00646         work2[iColumn]=w1;
00647         work3[iColumn]=work3Value;
00648         work4[iColumn]=work4Value;
00649       } else {
00650         work1[iColumn]=0.0;
00651         work2[iColumn]=0.0;
00652         work3[iColumn]=0.0;
00653         work4[iColumn]=0.0;
00654       } 
00655     } 
00656   } else {
00657     //iColumnust primal dual
00658     for (int iColumn=0;iColumn<numberTotal;iColumn++) {
00659       if (!flagged(iColumn)) {
00660         double z1=-zVec[iColumn];
00661         double w1=-wVec[iColumn];
00662         double work3Value=0.0;
00663         double work4Value=0.0;
00664         double directionElement=weights[iColumn];
00665         double value=primal[iColumn];
00666         if (directionNorm<fabs(directionElement)) {
00667           directionNorm=fabs(directionElement);
00668         } 
00669         //below does not feel right - can it be simplified because
00670         //of zero values for zVec and wVec
00671         if (lowerBound(iColumn)||
00672                         upperBound(iColumn)) {
00673           if (lowerBound(iColumn)) {
00674             double gap=lowerSlack[iColumn]+extra;
00675             double delta;
00676             if (!fakeLower(iColumn)) {
00677               delta = lower[iColumn]+lowerSlack[iColumn]
00678              -  value - directionElement;
00679             } else {
00680               delta = - directionElement;
00681             } 
00682             z1+=(mu_+zVec[iColumn]*delta)/gap;
00683             work3Value=-delta;
00684             if (lowerSlack[iColumn]<maximumPrimalStep*delta) {
00685               maximumPrimalStep=lowerSlack[iColumn]/delta;
00686               chosenPrimalSequence=iColumn;
00687             } 
00688             if (zVec[iColumn]>tolerance) {
00689               if (zVec[iColumn]<-z1*maximumDualStep) {
00690                 maximumDualStep=-zVec[iColumn]/z1;
00691                 chosenDualSequence=iColumn;
00692               } 
00693             } 
00694           } 
00695           if (upperBound(iColumn)) {
00696             double gap=upperSlack[iColumn]+extra;
00697             double delta;
00698             if (!fakeUpper(iColumn)) {
00699               delta = -upper[iColumn]+upperSlack[iColumn]
00700                         + value + directionElement;
00701             } else {
00702               delta = directionElement;
00703             } 
00704             w1+=(mu_+wVec[iColumn]*delta)/gap;
00705             work4Value=-delta;
00706             if (upperSlack[iColumn]<maximumPrimalStep*delta) {
00707               maximumPrimalStep=upperSlack[iColumn]/delta;
00708               chosenPrimalSequence=iColumn;
00709             } 
00710             if (wVec[iColumn]>tolerance) {
00711               if (wVec[iColumn]<-w1*maximumDualStep) {
00712                 maximumDualStep=-wVec[iColumn]/w1;
00713                 chosenDualSequence=iColumn;
00714               } 
00715             } 
00716           } 
00717         } else {
00718           //free
00719           double gap=fabs(value);
00720           double multiplier=1.0/gap;
00721           if (gap<1.0) {
00722             multiplier=1,0;
00723           } 
00724           z1+=multiplier*(mu_-directionElement*zVec[iColumn]);
00725           w1+=multiplier*(mu_+directionElement*wVec[iColumn]);
00726         } 
00727         work1[iColumn]=z1;
00728         work2[iColumn]=w1;
00729         work3[iColumn]=work3Value;
00730         work4[iColumn]=work4Value;
00731       } else {
00732         work1[iColumn]=0.0;
00733         work2[iColumn]=0.0;
00734         work3[iColumn]=0.0;
00735         work4[iColumn]=0.0;
00736       } 
00737     } 
00738   } 
00739   actualPrimalStep_=stepLength_*maximumPrimalStep;
00740   if (phase>0&&actualPrimalStep_>1.0) {
00741     actualPrimalStep_=1.0;
00742   } 
00743   actualDualStep_=stepLength_*maximumDualStep;
00744   if (phase>0&&actualDualStep_>1.0) {
00745     actualDualStep_=1.0;
00746   } 
00747   return directionNorm;
00748 }
00749 // findDirectionVector.
00750 double ClpPredictorCorrector::findDirectionVector(const int phase)
00751 {
00752   double projectionTolerance=projectionTolerance_;
00753   //temporary
00754   //projectionTolerance=1.0e-15;
00755   double errorCheck=0.9*maximumRHSError_/solutionNorm_;
00756   if (errorCheck>primalTolerance()) {
00757     if (errorCheck<projectionTolerance) {
00758       projectionTolerance=errorCheck;
00759     } 
00760   } else {
00761     if (primalTolerance()<projectionTolerance) {
00762       projectionTolerance=primalTolerance();
00763     } 
00764   } 
00765   double * newError = new double [numberRows_];
00766   double * work2 = deltaW_;
00767   int iColumn;
00768   int numberTotal = numberRows_+numberColumns_;
00769   //if flagged then entries zero so can do
00770   for (iColumn=0;iColumn<numberTotal;iColumn++)
00771     weights_[iColumn] = work2[iColumn] - solution_[iColumn];
00772   multiplyAdd(weights_+numberColumns_,numberRows_,-1.0,updateRegion_,0.0);
00773   matrix_->times(1.0,weights_,updateRegion_);
00774   bool goodSolve=false;
00775   double * regionSave=NULL;//for refinement
00776   int numberTries=0;
00777   double relativeError=COIN_DBL_MAX;
00778   double tryError=1.0e31;
00779   while (!goodSolve) {
00780     double lastError=relativeError;
00781     goodSolve=true;
00782     double maximumRHS = maximumAbsElement(updateRegion_,numberRows_);
00783     double scale=1.0;
00784     double unscale=1.0;
00785     if (maximumRHS>1.0e-30) {
00786       if (maximumRHS<=0.5) {
00787         double factor=2.0;
00788         while (maximumRHS<=0.5) {
00789           maximumRHS*=factor;
00790           scale*=factor;
00791         } /* endwhile */
00792       } else if (maximumRHS>=2.0) {
00793         double factor=0.5;
00794         while (maximumRHS>=2.0) {
00795           maximumRHS*=factor;
00796           scale*=factor;
00797         } /* endwhile */
00798       } 
00799       unscale=diagonalScaleFactor_/scale;
00800     } else {
00801       //effectively zero
00802       scale=0.0;
00803       unscale=0.0;
00804     } 
00805     //printf("--putting scales to 1.0\n");
00806     //scale=1.0;
00807     //unscale=1.0;
00808     multiplyAdd(NULL,numberRows_,0.0,updateRegion_,scale);
00809     cholesky_->solve(updateRegion_);
00810     multiplyAdd(NULL,numberRows_,0.0,updateRegion_,unscale);
00811     if (numberTries) {
00812       //refine?
00813       double scaleX=1.0;
00814       if (lastError>1.0e-5) 
00815         scaleX=0.8;
00816       multiplyAdd(regionSave,numberRows_,1.0,updateRegion_,scaleX);
00817     } 
00818     numberTries++;
00819     CoinZeroN(newError,numberRows_);
00820     multiplyAdd(updateRegion_,numberRows_,-1.0,weights_+numberColumns_,0.0);
00821     CoinZeroN(weights_,numberColumns_);
00822     matrix_->transposeTimes(1.0,updateRegion_,weights_);
00823     //if flagged then entries zero so can do
00824     for (iColumn=0;iColumn<numberTotal;iColumn++)
00825       weights_[iColumn] = weights_[iColumn]*diagonal_[iColumn]
00826         -work2[iColumn];
00827     multiplyAdd(weights_+numberColumns_,numberRows_,-1.0,newError,1.0);
00828     matrix_->times(1.0,weights_,newError);
00829     
00830     //now add in old Ax - doing extra checking
00831     double maximumRHSError=0.0;
00832     double maximumRHSChange=0.0;
00833     int iRow;
00834     char * dropped = cholesky_->rowsDropped();
00835     for (iRow=0;iRow<numberRows_;iRow++) {
00836       if (!dropped[iRow]) {
00837         double newValue=newError[iRow];
00838         double oldValue=errorRegion_[iRow];
00839         //severity of errors depend on signs
00840         //**later                                                             */
00841         if (fabs(newValue)>maximumRHSChange) {
00842           maximumRHSChange=fabs(newValue);
00843         } 
00844         double result=newValue+oldValue;
00845         if (fabs(result)>maximumRHSError) {
00846           maximumRHSError=fabs(result);
00847         } 
00848         newError[iRow]=result;
00849       } else {
00850         double newValue=newError[iRow];
00851         double oldValue=errorRegion_[iRow];
00852         if (fabs(newValue)>maximumRHSChange) {
00853           maximumRHSChange=fabs(newValue);
00854         } 
00855         double result=newValue+oldValue;
00856         newError[iRow]=result;
00857         //newError[iRow]=0.0;
00858         assert(updateRegion_[iRow]==0.0);
00859       } 
00860     } 
00861     relativeError = maximumRHSError/solutionNorm_;
00862     if (relativeError>tryError) 
00863       relativeError=tryError;
00864     if (relativeError<lastError) {
00865       maximumRHSChange_= maximumRHSChange;
00866       if (relativeError>1.0e-9
00867           ||numberTries>1) {
00868         handler_->message(CLP_BARRIER_ACCURACY,messages_)
00869           <<phase<<numberTries<<relativeError
00870           <<CoinMessageEol;
00871       } 
00872       if (relativeError>projectionTolerance&&numberTries<=3) {
00873         //try and refine
00874         goodSolve=false;
00875       } 
00876       //*** extra test here
00877       if (!goodSolve) {
00878         if (!regionSave) {
00879           regionSave = new double [numberRows_];
00880         } 
00881         CoinMemcpyN(updateRegion_,numberRows_,regionSave);
00882         multiplyAdd(newError,numberRows_,-1.0,updateRegion_,0.0);
00883       } 
00884     } else {
00885       //std::cout <<" worse residual = "<<relativeError;
00886       //bring back previous
00887       relativeError=lastError;
00888       CoinMemcpyN(regionSave,numberRows_,updateRegion_);
00889       multiplyAdd(updateRegion_,numberRows_,-1.0,weights_+numberColumns_,0.0);
00890       CoinZeroN(weights_,numberColumns_);
00891       matrix_->transposeTimes(1.0,updateRegion_,weights_);
00892       //if flagged then entries zero so can do
00893       for (iColumn=0;iColumn<numberTotal;iColumn++)
00894         weights_[iColumn] = weights_[iColumn]*diagonal_[iColumn]
00895           -work2[iColumn];
00896     } 
00897   } /* endwhile */
00898   delete [] regionSave;
00899   delete [] newError;
00900   return relativeError;
00901 }
00902 // createSolution.  Creates solution from scratch
00903 void ClpPredictorCorrector::createSolution()
00904 {
00905   int numberTotal = numberRows_+numberColumns_;
00906   int iColumn;
00907   double tolerance = primalTolerance();
00908   for (iColumn=0;iColumn<numberTotal;iColumn++) {
00909     if (upper_[iColumn]-lower_[iColumn]>tolerance) 
00910       clearFixed(iColumn);
00911     else
00912       setFixed(iColumn);
00913   }
00914   double initialValue =1.0e-12;
00915 
00916   double maximumObjective=1.0e-12;
00917   double objectiveNorm2=0.0;
00918   getNorms(cost_,numberTotal,maximumObjective,objectiveNorm2);
00919   if (maximumObjective<1.0e-12) {
00920     maximumObjective=1.0e-12;
00921   } 
00922   objectiveNorm2=sqrt(objectiveNorm2)/(double) numberTotal;
00923   objectiveNorm_=maximumObjective;
00924   scaleFactor_=1.0;
00925   if (maximumObjective>0.0) {
00926     if (maximumObjective<1.0) {
00927       scaleFactor_=maximumObjective;
00928     } else if (maximumObjective>1.0e4) {
00929       scaleFactor_=maximumObjective/1.0e4;
00930     } 
00931   } 
00932   if (scaleFactor_!=1.0) {
00933     objectiveNorm2*=scaleFactor_;
00934     multiplyAdd(NULL,numberTotal,0.0,cost_,1.0/scaleFactor_);
00935     objectiveNorm_=maximumObjective/scaleFactor_;
00936   } 
00937   baseObjectiveNorm_=objectiveNorm_;
00938   //accumulate fixed in dj region (as spare)
00939   //acumulate primal solution in primal region
00940   //DZ in lowerDual
00941   //DW in upperDual
00942   double infiniteCheck=1.0e40;
00943   //double     fakeCheck=1.0e10;
00944   //use weights region for work region
00945   for (iColumn=0;iColumn<numberTotal;iColumn++) {
00946     double primalValue = solution_[iColumn];
00947     clearFlagged(iColumn);
00948     clearFixedOrFree(iColumn);
00949     clearLowerBound(iColumn);
00950     clearUpperBound(iColumn);
00951     clearFakeLower(iColumn);
00952     clearFakeUpper(iColumn);
00953     if (!fixed(iColumn)) {
00954       dj_[iColumn]=0.0;
00955       diagonal_[iColumn]=1.0;
00956       weights_[iColumn]=1.0;
00957       double lowerValue=lower_[iColumn];
00958       double upperValue=upper_[iColumn];
00959 #if 0
00960       //fake it
00961       if (lowerValue<-fakeCheck&&upperValue>fakeCheck) {
00962         lowerValue=-1.0e5;
00963         upperValue=1.0e5;
00964         lower_[iColumn]=lowerValue;
00965         upper_[iColumn]=upperValue;
00966         std::cout<<"faking free variable "<<iColumn<<std::endl;
00967       } 
00968 #endif
00969       if (lowerValue>-infiniteCheck) {
00970         if (upperValue<infiniteCheck) {
00971           //upper and lower bounds
00972           setLowerBound(iColumn);
00973           setUpperBound(iColumn);
00974           if (lowerValue>=0.0) {
00975             solution_[iColumn]=lowerValue;
00976           } else if (upperValue<=0.0) {
00977             solution_[iColumn]=upperValue;
00978           } else {
00979             solution_[iColumn]=0.0;
00980           } 
00981         } else {
00982           //just lower bound
00983           setLowerBound(iColumn);
00984           if (lowerValue>=0.0) {
00985             solution_[iColumn]=lowerValue;
00986           } else {
00987             solution_[iColumn]=0.0;
00988           } 
00989         } 
00990       } else {
00991         if (upperValue<infiniteCheck) {
00992           //just upper bound
00993           setUpperBound(iColumn);
00994           if (upperValue<=0.0) {
00995             solution_[iColumn]=upperValue;
00996           } else {
00997             solution_[iColumn]=0.0;
00998           } 
00999         } else {
01000           //free
01001           setFixedOrFree(iColumn);
01002           solution_[iColumn]=0.0;
01003           //std::cout<<" free "<<i<<std::endl;
01004         } 
01005       } 
01006     } else {
01007       setFlagged(iColumn);
01008       setFixedOrFree(iColumn);
01009       setLowerBound(iColumn);
01010       setUpperBound(iColumn);
01011       dj_[iColumn]=primalValue;;
01012       solution_[iColumn]=lower_[iColumn];
01013       diagonal_[iColumn]=0.0;
01014       weights_[iColumn]=0.0;
01015     } 
01016   } 
01017   //   modify fixed RHS
01018   multiplyAdd(solution_+numberColumns_,numberRows_,1.0,errorRegion_,0.0);
01019   multiplyAdd(dj_+numberColumns_,numberRows_,-1.0,rhsFixRegion_,0.0);
01020   matrix_->times(-1.0,solution_,errorRegion_);
01021   //   create plausible RHS?
01022   matrix_->times(-1.0,dj_,rhsFixRegion_);
01023   rhsNorm_=maximumAbsElement(errorRegion_,numberRows_);
01024   if (rhsNorm_<1.0) {
01025     rhsNorm_=1.0;
01026   } 
01027   int * rowsDropped = new int [numberRows_];
01028   cholesky_->factorize(diagonal_,rowsDropped);
01029   if (cholesky_->status()) {
01030     std::cout << "singular on initial cholesky?" <<std::endl;
01031     cholesky_->resetRowsDropped();
01032     //cholesky_->factorize(rowDropped_);
01033     //if (cholesky_->status()) {
01034       //std::cout << "bad cholesky??? (after retry)" <<std::endl;
01035       //abort();
01036     //} 
01037   } 
01038   delete [] rowsDropped;
01039   cholesky_->solve(errorRegion_);
01040   //create information for solution
01041   multiplyAdd(errorRegion_,numberRows_,-1.0,weights_+numberColumns_,0.0);
01042   CoinZeroN(weights_,numberColumns_);
01043   matrix_->transposeTimes(1.0,errorRegion_,weights_);
01044   //do reduced costs
01045   CoinZeroN(dj_+numberColumns_,numberRows_);
01046   for ( iColumn=0;iColumn<numberColumns_;iColumn++) {
01047     double value=cost_[iColumn];
01048     if (lowerBound(iColumn)) {
01049       value+=linearPerturbation_;
01050     } 
01051     if (upperBound(iColumn)) {
01052       value-=linearPerturbation_;
01053     } 
01054     dj_[iColumn]=value;
01055   } 
01056   initialValue=1.0e2;
01057   if (rhsNorm_*1.0e-2>initialValue) {
01058     initialValue=rhsNorm_*1.0e-2;
01059   } 
01060   double smallestBoundDifference=COIN_DBL_MAX;
01061   double * fakeSolution = weights_;
01062   for ( iColumn=0;iColumn<numberTotal;iColumn++) {
01063     if (!flagged(iColumn)) {
01064       if (lower_[iColumn]-fakeSolution[iColumn]>initialValue) {
01065         initialValue=lower_[iColumn]-fakeSolution[iColumn];
01066       } 
01067       if (fakeSolution[iColumn]-upper_[iColumn]>initialValue) {
01068         initialValue=fakeSolution[iColumn]-upper_[iColumn];
01069       } 
01070       if (upper_[iColumn]-lower_[iColumn]<smallestBoundDifference) {
01071         smallestBoundDifference=upper_[iColumn]-lower_[iColumn];
01072       } 
01073     } 
01074   } 
01075   solutionNorm_=1.0e-12;
01076   handler_->message(CLP_BARRIER_SAFE,messages_)
01077     <<initialValue<<objectiveNorm_
01078     <<CoinMessageEol;
01079   long strategy=1;
01080   double extra=1.0e-10;
01081   double largeGap=1.0e15;
01082   double safeObjectiveValue=2.0*objectiveNorm_;
01083   double safeFree=1.0e-1*initialValue;
01084   safeFree=1.0;
01085   double zwLarge=1.0e2*initialValue;
01086   //zwLarge=1.0e40;
01087   for ( iColumn=0;iColumn<numberTotal;iColumn++) {
01088     if (!flagged(iColumn)) {
01089       double lowerValue=lower_[iColumn];
01090       double upperValue=upper_[iColumn];
01091       double objectiveValue = cost_[iColumn];
01092       double newValue;
01093       double low;
01094       double high;
01095       double randomZ =0.5*CoinDrand48()+0.5;
01096       double randomW =0.5*CoinDrand48()+0.5;
01097       double setPrimal=initialValue;
01098       if (strategy==0) {
01099         randomZ=1.0;
01100         randomW=1.0;
01101       } 
01102       if (lowerBound(iColumn)) {
01103         if (upperBound(iColumn)) {
01104           //upper and lower bounds
01105           if (upperValue-lowerValue>2.0*setPrimal) {
01106             double fakeValue=fakeSolution[iColumn];
01107             if (fakeValue<lowerValue+setPrimal) {
01108               fakeValue=lowerValue+setPrimal;
01109             } 
01110             if (fakeValue>upperValue-setPrimal) {
01111               fakeValue=upperValue-setPrimal;
01112             } 
01113             newValue=fakeValue;
01114             low=fakeValue-lowerValue;
01115             high=upperValue-fakeValue;
01116           } else {
01117             newValue=0.5*(upperValue+lowerValue);
01118             low=setPrimal;
01119             high=setPrimal;
01120           } 
01121           low *= randomZ;
01122           high *= randomW;
01123           double s = low+extra;
01124           double ratioZ;
01125           if (s<zwLarge) {
01126             ratioZ=1.0;
01127           } else {
01128             ratioZ=sqrt(zwLarge/s);
01129           } 
01130           double t = high+extra;
01131           double ratioW;
01132           if (t<zwLarge) {
01133             ratioW=1.0;
01134           } else {
01135             ratioW=sqrt(zwLarge/t);
01136           } 
01137           //modify s and t
01138           if (s>largeGap) {
01139             s=largeGap;
01140           } 
01141           if (t>largeGap) {
01142             t=largeGap;
01143           } 
01144           //modify if long long way away from bound
01145           if (objectiveValue>=0.0) {
01146             zVec_[iColumn]=objectiveValue + safeObjectiveValue*ratioZ;
01147             wVec_[iColumn]=safeObjectiveValue*ratioW;
01148           } else {
01149             zVec_[iColumn]=safeObjectiveValue*ratioZ;
01150             wVec_[iColumn]=-objectiveValue + safeObjectiveValue*ratioW;
01151           } 
01152           diagonal_[iColumn] = (t*s)/(s*wVec_[iColumn]+t*zVec_[iColumn]);
01153         } else {
01154           //just lower bound
01155           double fakeValue=fakeSolution[iColumn];
01156           if (fakeValue<lowerValue+setPrimal) {
01157             fakeValue=lowerValue+setPrimal;
01158           } 
01159           newValue=fakeValue;
01160           low=fakeValue-lowerValue;
01161           low *= randomZ;
01162           high=0.0;
01163           double s = low+extra;
01164           double ratioZ;
01165           if (s<zwLarge) {
01166             ratioZ=1.0;
01167           } else {
01168             ratioZ=sqrt(zwLarge/s);
01169           } 
01170           //modify s
01171           if (s>largeGap) {
01172             s=largeGap;
01173           } 
01174           if (objectiveValue>=0.0) {
01175             zVec_[iColumn]=objectiveValue + safeObjectiveValue*ratioZ;
01176             wVec_[iColumn]=0.0;
01177           } else {
01178             zVec_[iColumn]=safeObjectiveValue*ratioZ;
01179             wVec_[iColumn]=0.0;
01180           } 
01181           diagonal_[iColumn] = s/zVec_[iColumn];
01182         } 
01183       } else {
01184         if (upperBound(iColumn)) {
01185           //just upper bound
01186           double fakeValue=fakeSolution[iColumn];
01187           if (fakeValue>upperValue-setPrimal) {
01188             fakeValue=upperValue-setPrimal;
01189           } 
01190           newValue=fakeValue;
01191           low=0.0;
01192           high=upperValue-fakeValue;
01193           high*=randomW;
01194           double t = high+extra;
01195           double ratioW;
01196           if (t<zwLarge) {
01197             ratioW=1.0;
01198           } else {
01199             ratioW=sqrt(zwLarge/t);
01200           } 
01201           //modify t
01202           if (t>largeGap) {
01203             t=largeGap;
01204           } 
01205           if (objectiveValue>=0.0) {
01206             zVec_[iColumn]=0.0;
01207             wVec_[iColumn]=safeObjectiveValue*ratioW;
01208           } else {
01209             zVec_[iColumn]=0.0;
01210             wVec_[iColumn]=-objectiveValue + safeObjectiveValue*ratioW;
01211           } 
01212           diagonal_[iColumn] =  t/wVec_[iColumn];
01213         } else {
01214           //free
01215           zVec_[iColumn]=safeObjectiveValue;
01216           wVec_[iColumn]=safeObjectiveValue;
01217           newValue=fakeSolution[iColumn];
01218           if (newValue>=0.0) {
01219             if (newValue<safeFree) {
01220               newValue=safeFree;
01221             } 
01222           } else {
01223             if (newValue>-safeFree) {
01224               newValue=-safeFree;
01225             } 
01226           } 
01227           if (fabs(newValue)>1.0) {
01228             diagonal_[iColumn]=fabs(newValue)/(zVec_[iColumn]+wVec_[iColumn]);
01229           } else {
01230             diagonal_[iColumn]=1.0/(zVec_[iColumn]+wVec_[iColumn]);
01231           } 
01232           low=0.0;
01233           high=0.0;
01234         } 
01235       } 
01236       lowerSlack_[iColumn]=low;
01237       upperSlack_[iColumn]=high;
01238       solution_[iColumn]=newValue;
01239     } else {
01240       lowerSlack_[iColumn]=0.0;
01241       upperSlack_[iColumn]=0.0;
01242       solution_[iColumn]=lower_[iColumn];
01243       zVec_[iColumn]=0.0;
01244       wVec_[iColumn]=0.0;
01245       diagonal_[iColumn]=0.0;
01246     } 
01247   } 
01248   solutionNorm_ =  maximumAbsElement(solution_,numberTotal);
01249 }
01250 // complementarityGap.  Computes gap
01251 //phase 0=as is , 1 = after predictor , 2 after corrector
01252 double ClpPredictorCorrector::complementarityGap(int & numberComplementarityPairs,
01253                           const int phase)
01254 {
01255   double gap=0.0;
01256   //seems to be same coding for phase = 1 or 2
01257   numberComplementarityPairs=0;
01258   double toleranceGap=0.0;
01259   double largestGap=0.0;
01260   //seems to be same coding for phase = 1 or 2
01261   int numberNegativeGaps=0;
01262   double sumNegativeGap=0.0;
01263   double largeGap=1.0e2*solutionNorm_;
01264   if (largeGap<1.0e10) {
01265     largeGap=1.0e10;
01266   } 
01267   largeGap=1.0e30;
01268   double dualTolerance =  dblParam_[ClpDualTolerance];
01269   double primalTolerance =  dblParam_[ClpPrimalTolerance];
01270   dualTolerance=dualTolerance/scaleFactor_;
01271   double * zVec = zVec_;
01272   double * wVec = wVec_;
01273   double * primal = solution_;
01274   double * lower = lower_;
01275   double * upper = upper_;
01276   double * lowerSlack = lowerSlack_;
01277   double * upperSlack = upperSlack_;
01278   double * work1 = deltaZ_;
01279   double * work2 = deltaW_;
01280   double * weights = weights_;
01281   for (int iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
01282     if (!fixedOrFree(iColumn)) {
01283       //can collapse as if no lower bound both zVec and work1 0.0
01284       if (lowerBound(iColumn)) {
01285         double dualValue;
01286         double primalValue;
01287         if (!phase) {
01288           dualValue=zVec[iColumn];
01289           primalValue=lowerSlack[iColumn];
01290         } else {
01291           double change;
01292           if (!fakeLower(iColumn)) {
01293             change =primal[iColumn]+weights[iColumn]-lowerSlack[iColumn]-lower[iColumn];
01294           } else {
01295             change =weights[iColumn];
01296           } 
01297           dualValue=zVec[iColumn]+actualDualStep_*work1[iColumn];
01298           primalValue=lowerSlack[iColumn]+actualPrimalStep_*change;
01299         } 
01300         //reduce primalValue
01301         if (primalValue>largeGap) {
01302           primalValue=largeGap;
01303         } 
01304         double gapProduct=dualValue*primalValue;
01305         if (gapProduct<0.0) {
01306           //cout<<"negative gap component "<<iColumn<<" "<<dualValue<<" "<<
01307               //primalValue<<endl;
01308           numberNegativeGaps++;
01309           sumNegativeGap-=gapProduct;
01310           gapProduct=0.0;
01311         } 
01312         gap+=gapProduct;
01313         if (gapProduct>largestGap) {
01314           largestGap=gapProduct;
01315         } 
01316         if (dualValue>dualTolerance&&primalValue>primalTolerance) {
01317           toleranceGap+=dualValue*primalValue;
01318         } 
01319         numberComplementarityPairs++;
01320       } 
01321       if (upperBound(iColumn)) {
01322         double dualValue;
01323         double primalValue;
01324         if (!phase) {
01325           dualValue=wVec[iColumn];
01326           primalValue=upperSlack[iColumn];
01327         } else {
01328           double change;
01329           if (!fakeUpper(iColumn)) {
01330             change =upper[iColumn]-primal[iColumn]-weights[iColumn]-upperSlack[iColumn];
01331           } else {
01332             change =-weights[iColumn];
01333           } 
01334           dualValue=wVec[iColumn]+actualDualStep_*work2[iColumn];
01335           primalValue=upperSlack[iColumn]+actualPrimalStep_*change;
01336         } 
01337         //reduce primalValue
01338         if (primalValue>largeGap) {
01339           primalValue=largeGap;
01340         } 
01341         double gapProduct=dualValue*primalValue;
01342         if (gapProduct<0.0) {
01343           //cout<<"negative gap component "<<iColumn<<" "<<dualValue<<" "<<
01344               //primalValue<<endl;
01345           numberNegativeGaps++;
01346           sumNegativeGap-=gapProduct;
01347           gapProduct=0.0;
01348         } 
01349         gap+=gapProduct;
01350         if (gapProduct>largestGap) {
01351           largestGap=gapProduct;
01352         } 
01353         if (dualValue>dualTolerance&&primalValue>primalTolerance) {
01354           toleranceGap+=dualValue*primalValue;
01355         } 
01356         numberComplementarityPairs++;
01357       } 
01358     } 
01359   } 
01360   if (!phase&&numberNegativeGaps) {
01361       handler_->message(CLP_BARRIER_NEGATIVE_GAPS,messages_)
01362     <<numberNegativeGaps<<sumNegativeGap
01363     <<CoinMessageEol;
01364   } 
01365   //in case all free!
01366   if (!numberComplementarityPairs) {
01367     numberComplementarityPairs=1;
01368   } 
01369   return gap;
01370 }
01371 // setupForSolve.
01372 //phase 0=affine , 1 = corrector , 2 = primal-dual
01373 void ClpPredictorCorrector::setupForSolve(const int phase)
01374 {
01375   double extra =eExtra;
01376   double * zVec = zVec_;
01377   double * wVec = wVec_;
01378   double * primal = solution_;
01379   double * dual = dj_;
01380   double * lower = lower_;
01381   double * upper = upper_;
01382   double * lowerSlack = lowerSlack_;
01383   double * upperSlack = upperSlack_;
01384   double * work2 = deltaW_;
01385   double * work3 = deltaS_;
01386   double * work4 = deltaT_;
01387   double * diagonal = diagonal_;
01388 
01389   int iColumn;
01390   switch (phase) {
01391   case 0:
01392     for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
01393       if (!flagged(iColumn)) {
01394         double value= dual[iColumn];
01395         if (lowerBound(iColumn)) {
01396           if (!fakeLower(iColumn)) {
01397             value+=zVec[iColumn]*
01398               (primal[iColumn]-lowerSlack[iColumn]-lower[iColumn])/
01399                  (lowerSlack[iColumn]+extra);
01400           } 
01401         } 
01402         if (upperBound(iColumn)) {
01403           if (!fakeUpper(iColumn)) {
01404             value+=wVec[iColumn]*
01405               (primal[iColumn]+upperSlack[iColumn]-upper[iColumn])/
01406                  (upperSlack[iColumn]+extra);
01407           } 
01408         } 
01409         work2[iColumn]=diagonal[iColumn]*value;
01410       } else {
01411         work2[iColumn]=0.0;
01412       } 
01413     } 
01414     break;
01415   case 1:
01416     for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
01417       if (!flagged(iColumn)) {
01418         double value= 0.0;
01419         if (lowerBound(iColumn)) {
01420           if (!fakeLower(iColumn)) {
01421             value-=(mu_-work3[iColumn]-zVec[iColumn]*
01422               (primal[iColumn]-lowerSlack[iColumn]-lower[iColumn]))/
01423                  (lowerSlack[iColumn]+extra);
01424           } else {
01425             value-=(mu_-work3[iColumn])/(lowerSlack[iColumn]+extra);
01426           } 
01427         } 
01428         if (upperBound(iColumn)) {
01429           if (!fakeUpper(iColumn)) {
01430             value+=(mu_-work4[iColumn]+wVec[iColumn]*
01431               (primal[iColumn]+upperSlack[iColumn]-upper[iColumn]))/
01432                  (upperSlack[iColumn]+extra);
01433           } else {
01434             value+=(mu_-work4[iColumn])/(upperSlack[iColumn]+extra);
01435           } 
01436         } 
01437         work2[iColumn]=diagonal[iColumn]*(dual[iColumn]+value);
01438       } else {
01439         work2[iColumn]=0.0;
01440       } 
01441     } 
01442     break;
01443   case 2:
01444     for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
01445       if (!flagged(iColumn)) {
01446         double value= 0.0;
01447         if (lowerBound(iColumn)) {
01448           if (!fakeLower(iColumn)) {
01449             value-=(mu_-zVec[iColumn]*
01450               (primal[iColumn]-lowerSlack[iColumn]-lower[iColumn]))/
01451                  (lowerSlack[iColumn]+extra);
01452           } else {
01453             value-=(mu_-zVec[iColumn])/ (lowerSlack[iColumn]+extra);
01454           } 
01455         } 
01456         if (upperBound(iColumn)) {
01457           if (!fakeUpper(iColumn)) {
01458             value+=(mu_+wVec[iColumn]*
01459               (primal[iColumn]+upperSlack[iColumn]-upper[iColumn]))/
01460                  (upperSlack[iColumn]+extra);
01461           } else {
01462             value+=(mu_+wVec[iColumn])/ (upperSlack[iColumn]+extra);
01463           } 
01464         } 
01465         work2[iColumn]=diagonal[iColumn]*(dual[iColumn]+value);
01466       } else {
01467         work2[iColumn]=0.0;
01468       } 
01469     } 
01470     break;
01471   } /* endswitch */
01472 }
01473 //method: sees if looks plausible change in complementarity
01474 bool ClpPredictorCorrector::checkGoodMove(const bool doCorrector)
01475 {
01476   const double beta3 = 0.99997;
01477   bool goodMove=false;
01478   int nextNumber;
01479   int numberTotal = numberRows_+numberColumns_;
01480   double nextGap=complementarityGap(nextNumber,2);
01481   double step;
01482   if (actualDualStep_>actualPrimalStep_) {
01483     step=actualDualStep_;
01484   } else {
01485     step=actualPrimalStep_;
01486   } 
01487   double testValue=1.0-step*(1.0-beta3);
01488   //testValue=0.0;
01489   testValue*=complementarityGap_;
01490   if (nextGap<testValue) {
01491     //std::cout <<"predicted duality gap "<<nextGap<<std::endl;
01492     goodMove=true;
01493   } else if(doCorrector) {
01494     //if (actualDualStep_<actualPrimalStep_) {
01495       //step=actualDualStep_;
01496     //} else {
01497       //step=actualPrimalStep_;
01498     //} 
01499     goodMove=checkGoodMove2(step);
01500   } else {
01501     goodMove=true;
01502   } 
01503   if (!goodMove) {
01504     //try smaller of two
01505     if (actualDualStep_<actualPrimalStep_) {
01506       step=actualDualStep_;
01507     } else {
01508       step=actualPrimalStep_;
01509     } 
01510     if (step>1.0) {
01511       step=1.0;
01512     } 
01513     actualPrimalStep_=step;
01514     actualDualStep_=step;
01515     while (!goodMove) {
01516       goodMove=checkGoodMove2(step);
01517       if (step<1.0e-10) {
01518         break;
01519       } 
01520       step*=0.5;
01521       actualPrimalStep_=step;
01522       actualDualStep_=step;
01523     } /* endwhile */
01524     if (doCorrector) {
01525       //say bad move if both small
01526       if (numberIterations_&1) {
01527         if (actualPrimalStep_<1.0e-2&&actualDualStep_<1.0e-2) {
01528           goodMove=false;
01529         } 
01530       } else {
01531         if (actualPrimalStep_<1.0e-5&&actualDualStep_<1.0e-5) {
01532           goodMove=false;
01533         } 
01534         if (actualPrimalStep_<1.0e-5*actualDualStep_<1.0e-20) {
01535           goodMove=false;
01536         } 
01537       } 
01538     } 
01539   } 
01540   if (goodMove) {
01541     //compute delta in objectives
01542     double deltaObjectivePrimal=0.0;
01543     double deltaObjectiveDual=
01544       innerProduct(updateRegion_,numberRows_,
01545                    rhsFixRegion_);
01546     double error=0.0;
01547     double * work3 = deltaS_;
01548     CoinZeroN(work3,numberColumns_);
01549     CoinMemcpyN(updateRegion_,numberRows_,work3+numberColumns_);
01550     matrix_->transposeTimes(-1.0,updateRegion_,work3);
01551     double * work1 = deltaZ_;
01552     double * work2 = deltaW_;
01553     double * lower = lower_;
01554     double * upper = upper_;
01555     //direction vector in weights
01556     double * weights = weights_;
01557     double * cost = cost_;
01558     //double sumPerturbCost=0.0;
01559     for (int iColumn=0;iColumn<numberTotal;iColumn++) {
01560       if (!flagged(iColumn)) {
01561         if (lowerBound(iColumn)) {
01562           //sumPerturbCost+=weights[iColumn];
01563           deltaObjectiveDual+=work1[iColumn]*lower[iColumn];
01564         } 
01565         if (upperBound(iColumn)) {
01566           //sumPerturbCost-=weights[iColumn];
01567           deltaObjectiveDual-=work2[iColumn]*upper[iColumn];
01568         } 
01569         double change = fabs(work3[iColumn]-work1[iColumn]+work2[iColumn]);
01570         error = max (change,error);
01571       } 
01572       deltaObjectivePrimal += cost[iColumn] * weights[iColumn];
01573     } 
01574     //deltaObjectivePrimal+=sumPerturbCost*linearPerturbation_;
01575     double testValue;
01576     if (error>0.0) {
01577       testValue=1.0e1*maximumDualError_/error;
01578     } else {
01579       testValue=1.0e1;
01580     } 
01581     if (testValue<actualDualStep_) {
01582       handler_->message(CLP_BARRIER_REDUCING,messages_)
01583       <<"dual"<<actualDualStep_
01584       << testValue
01585       <<CoinMessageEol;
01586       actualDualStep_=testValue;
01587     } 
01588   } 
01589   if (maximumRHSError_<1.0e1*solutionNorm_*primalTolerance()
01590                             &&maximumRHSChange_>1.0e-16*solutionNorm_) {
01591     //check change in AX not too much
01592     //??? could be dropped row going infeasible
01593     double ratio = 1.0e1*maximumRHSError_/maximumRHSChange_;
01594     if (ratio<actualPrimalStep_) {
01595       handler_->message(CLP_BARRIER_REDUCING,messages_)
01596       <<"primal"<<actualPrimalStep_
01597       <<ratio
01598       <<CoinMessageEol;
01599       if (ratio>1.0e-6) {
01600         actualPrimalStep_=ratio;
01601       } else {
01602         actualPrimalStep_=ratio;
01603         //std::cout <<"sign we should be stopping"<<std::endl;
01604       } 
01605     } 
01606   } 
01607   return goodMove;
01608 }
01609 //:  checks for one step size
01610 bool ClpPredictorCorrector::checkGoodMove2(const double move)
01611 {
01612   double complementarityMultiplier =1.0/numberComplementarityPairs_;
01613   const double gamma = 1.0e-8;
01614   const double gammap = 1.0e-8;
01615   const double gammad = 1.0e-8;
01616   int nextNumber;
01617   double nextGap=complementarityGap(nextNumber,2);
01618   double lowerBoundGap = gamma*nextGap*complementarityMultiplier;
01619   bool goodMove=true;
01620   double * deltaZ = deltaZ_;
01621   double * deltaW = deltaW_;
01622   double * deltaS = deltaS_;
01623   double * deltaT = deltaT_;
01624   double * zVec = zVec_;
01625   double * wVec = wVec_;
01626   double * lowerSlack = lowerSlack_;
01627   double * upperSlack = upperSlack_;
01628   for (int iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
01629     if (!flagged(iColumn)) {
01630       if (lowerBound(iColumn)) {
01631         double part1=lowerSlack[iColumn]+actualPrimalStep_*deltaS[iColumn];
01632         double part2=zVec[iColumn]+actualDualStep_*deltaZ[iColumn];
01633         if (part1*part2<lowerBoundGap) {
01634           goodMove=false;
01635           break;
01636         } 
01637       } 
01638       if (upperBound(iColumn)) {
01639         double part1=upperSlack[iColumn]+actualPrimalStep_*deltaT[iColumn];
01640         double part2=wVec[iColumn]+actualDualStep_*deltaW[iColumn];
01641         if (part1*part2<lowerBoundGap) {
01642           goodMove=false;
01643           break;
01644         } 
01645       } 
01646     } 
01647   } 
01648   //      Satisfy g_p(alpha)?
01649   if (rhsNorm_>solutionNorm_) {
01650     solutionNorm_=rhsNorm_;
01651   } 
01652   double errorCheck=maximumRHSError_/solutionNorm_;
01653   if (errorCheck<maximumBoundInfeasibility_) {
01654     errorCheck=maximumBoundInfeasibility_;
01655   } 
01656   //scale
01657   if ((1.0-move)*errorCheck>primalTolerance()) {
01658     if (nextGap<gammap*(1.0-move)*errorCheck) {
01659       goodMove=false;
01660     } 
01661   } 
01662   //      Satisfy g_d(alpha)?
01663   errorCheck=maximumDualError_/objectiveNorm_;
01664   if ((1.0-move)*errorCheck>dualTolerance()) {
01665     if (nextGap<gammad*(1.0-move)*errorCheck) {
01666       goodMove=false;
01667     } 
01668   } 
01669   return goodMove;
01670 }
01671 // updateSolution.  Updates solution at end of iteration
01672 //returns number fixed
01673 int ClpPredictorCorrector::updateSolution()
01674 {
01675   //update pi
01676   multiplyAdd(updateRegion_,numberRows_,actualDualStep_,dual_,1.0);
01677   CoinZeroN(errorRegion_,numberRows_);
01678   CoinZeroN(rhsFixRegion_,numberRows_);
01679   double maximumRhsInfeasibility=0.0;
01680   double maximumBoundInfeasibility=0.0;
01681   double maximumDualError=1.0e-12;
01682   double primalObjectiveValue=0.0;
01683   double dualObjectiveValue=0.0;
01684   double solutionNorm=1.0e-12;
01685   int numberKilled=0;
01686   double freeMultiplier=1.0e6;
01687   double trueNorm =diagonalNorm_/diagonalScaleFactor_;
01688   if (freeMultiplier<trueNorm) {
01689     freeMultiplier=trueNorm;
01690   } 
01691   if (freeMultiplier>1.0e12) {
01692     freeMultiplier=1.0e12;
01693   } 
01694   freeMultiplier=0.5/freeMultiplier;
01695   double condition = cholesky_->choleskyCondition();
01696   bool caution;
01697   if ((condition<1.0e10&&trueNorm<1.0e12)||numberIterations_<20) {
01698     caution=false;
01699   } else {
01700     caution=true;
01701   } 
01702   // do reduced costs
01703   CoinMemcpyN(dual_,numberRows_,dj_+numberColumns_);
01704   CoinMemcpyN(cost_,numberColumns_,dj_);
01705   matrix_->transposeTimes(-1.0,dual_,dj_);
01706   double extra=eExtra;
01707   double largeGap=1.0e2*solutionNorm_;
01708   if (largeGap<1.0e2) {
01709     largeGap=1.0e2;
01710   } 
01711   double dualFake=0.0;
01712   double dualTolerance =  dblParam_[ClpDualTolerance];
01713   dualTolerance=dualTolerance/scaleFactor_;
01714   if (dualTolerance<1.0e-12) {
01715     dualTolerance=1.0e-12;
01716   } 
01717   double offsetObjective=0.0;
01718   const double killTolerance=primalTolerance();
01719   double qDiagonal;
01720   if (mu_<1.0) {
01721     qDiagonal=1.0e-8;
01722   } else {
01723     qDiagonal=1.0e-8*mu_;
01724   } 
01725   double norm=1.0e-12;
01726   double widenGap=1.0e1;
01727   //largest allowable ratio of lowerSlack/zVec (etc)
01728   double largestRatio;
01729   double epsilonBase;
01730   double diagonalLimit;
01731   if (!caution) {
01732     epsilonBase=eBase;
01733     largestRatio=eRatio;
01734     diagonalLimit=eDiagonal;
01735   } else {
01736     epsilonBase=eBaseCaution;
01737     largestRatio=eRatioCaution;
01738     diagonalLimit=eDiagonalCaution;
01739   } 
01740   double smallGap=1.0e2;
01741   double maximumDJInfeasibility=0.0;
01742   int numberIncreased=0;
01743   int numberDecreased=0;
01744   double largestDiagonal=0.0;
01745   double smallestDiagonal=1.0e50;
01746   double * zVec = zVec_;
01747   double * wVec = wVec_;
01748   double * primal = solution_;
01749   double * dual = dj_;
01750   double * lower = lower_;
01751   double * upper = upper_;
01752   double * lowerSlack = lowerSlack_;
01753   double * upperSlack = upperSlack_;
01754   double * work1 = deltaZ_;
01755   double * work2 = deltaW_;
01756   double * diagonal = diagonal_;
01757   //direction vector in weights
01758   double * weights = weights_;
01759   double * cost = cost_;
01760   for (int iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
01761     if (!flagged(iColumn)) {
01762       double reducedCost=dual[iColumn];
01763       bool thisKilled=false;
01764       double zValue = zVec[iColumn] + actualDualStep_*work1[iColumn];
01765       double wValue = wVec[iColumn] + actualDualStep_*work2[iColumn];
01766       zVec[iColumn]=zValue;
01767       wVec[iColumn]=wValue;
01768       double thisWeight=weights[iColumn];
01769       double oldPrimal = primal[iColumn];
01770       double newPrimal=primal[iColumn]+actualPrimalStep_*thisWeight;
01771       double primalObjectiveThis=newPrimal*cost[iColumn];
01772       double dualObjectiveThis=0.0;
01773       double t=extra;
01774       double s=extra;
01775       double kill;
01776       if (fabs(newPrimal)>1.0e4) {
01777         kill=killTolerance*1.0e-4*newPrimal;
01778         //kill=killTolerance;
01779         //kill*=1.0e-3;//???????
01780       } else {
01781         kill=killTolerance;
01782         //kill*=1.0e-3;//???????
01783       } 
01784       kill*=1.0e-3;//be conservative
01785       double smallerSlack=COIN_DBL_MAX;
01786       double largerzw;
01787       if (zValue>wValue) {
01788         largerzw=zValue;
01789       } else {
01790         largerzw=wValue;
01791       } 
01792       if (lowerBound(iColumn)) {
01793         double oldSlack = lowerSlack[iColumn];
01794         double newSlack;
01795         if (!fakeLower(iColumn)) {
01796           newSlack=
01797                    lowerSlack[iColumn]+actualPrimalStep_*(oldPrimal-oldSlack
01798                 + thisWeight-lower[iColumn]);
01799         } else {
01800           newSlack= lowerSlack[iColumn]+actualPrimalStep_*thisWeight;
01801           if (newSlack<0.0) {
01802             abort();
01803           } 
01804         } 
01805         double epsilon = fabs(newSlack)*epsilonBase;
01806         if (epsilon>1.0e-5) {
01807           //cout<<"bad"<<endl;
01808           epsilon=1.0e-5;
01809         } 
01810         //for changing slack
01811         double zValue2 = zValue;
01812         if (zValue2<epsilonBase) {
01813           zValue2=epsilonBase;
01814         } 
01815         //make sure reasonable
01816         if (zValue<epsilon) {
01817           zValue=epsilon;
01818         } 
01819         //store modified zVec
01820         //no zVec[iColumn]=zValue;
01821         double feasibleSlack=newPrimal-lower[iColumn];
01822         if (feasibleSlack>0.0&&newSlack>0.0) {
01823           double smallGap2=smallGap;
01824           if (fabs(0.1*newPrimal)>smallGap) {
01825             smallGap2=0.1*fabs(newPrimal);
01826           } 
01827           if (!fakeLower(iColumn)) {
01828             double larger;
01829             if (newSlack>feasibleSlack) {
01830               larger=newSlack;
01831             } else {
01832               larger=feasibleSlack;
01833             } 
01834             if (fabs(feasibleSlack-newSlack)<1.0e-6*larger) {
01835               newSlack=feasibleSlack;
01836             } 
01837             //set FAKE here
01838             if (newSlack>zValue2*largestRatio&&newSlack>smallGap2) {
01839               setFakeLower(iColumn);
01840               newSlack=zValue2*largestRatio;
01841               if (newSlack<smallGap2) {
01842                 newSlack=smallGap2;
01843               } 
01844               numberDecreased++;
01845             } 
01846           } else {
01847             newSlack=zValue2*largestRatio;
01848             if (newSlack<smallGap2) {
01849               newSlack=smallGap2;
01850             } 
01851             if (newSlack>largeGap) {
01852               //increase up to smaller of z.. and largeGap
01853               newSlack=largeGap;
01854             } 
01855             if (newSlack>widenGap*oldSlack) {
01856               newSlack=widenGap*oldSlack;
01857               numberIncreased++;
01858               //cout<<"wider "<<j<<" "<<newSlack<<" "<<oldSlack<<" "<<feasibleSlack<<endl;
01859             } 
01860             if (newSlack>feasibleSlack) {
01861               newSlack=feasibleSlack;
01862               clearFakeLower(iColumn);
01863             } 
01864           } 
01865         } 
01866         if (zVec[iColumn]>dualTolerance) {
01867           dualObjectiveThis+=lower[iColumn]*zVec[iColumn];
01868         } 
01869         lowerSlack[iColumn]=newSlack;
01870         if (newSlack<smallerSlack) {
01871           smallerSlack=newSlack;
01872         } 
01873         if (!fakeLower(iColumn)) {
01874           double infeasibility = fabs(newPrimal-lowerSlack[iColumn]-lower[iColumn]);
01875           if (infeasibility>maximumBoundInfeasibility) {
01876             maximumBoundInfeasibility=infeasibility;
01877           } 
01878         } 
01879         if (lowerSlack[iColumn]<=kill&&fabs(newPrimal-lower[iColumn])<=kill) {
01880           //may be better to leave at value?
01881           newPrimal=lower[iColumn];
01882           lowerSlack[iColumn]=0.0;
01883           thisKilled=true;
01884           //cout<<j<<" l "<<reducedCost<<" "<<zVec[iColumn]<<endl;
01885         } else {
01886           s+=lowerSlack[iColumn];
01887         } 
01888       } 
01889       if (upperBound(iColumn)) {
01890         //reducedCost-=perturbation;
01891         //primalObjectiveThis-=perturbation*newPrimal;
01892         double oldSlack = upperSlack[iColumn];
01893         double newSlack;
01894         if (!fakeUpper(iColumn)) {
01895           newSlack=
01896                    upperSlack[iColumn]+actualPrimalStep_*(-oldPrimal-oldSlack
01897                 - thisWeight+upper[iColumn]);
01898         } else {
01899           newSlack= upperSlack[iColumn]-actualPrimalStep_*thisWeight;
01900         } 
01901         double epsilon = fabs(newSlack)*epsilonBase;
01902         if (epsilon>1.0e-5) {
01903           //cout<<"bad"<<endl;
01904           epsilon=1.0e-5;
01905         } 
01906         //for changing slack
01907         double wValue2 = wValue;
01908         if (wValue2<epsilonBase) {
01909           wValue2=epsilonBase;
01910         } 
01911         //make sure reasonable
01912         if (wValue<epsilon) {
01913           wValue=epsilon;
01914         } 
01915         //store modified wVec
01916         //no wVec[iColumn]=wValue;
01917         double feasibleSlack=upper[iColumn]-newPrimal;
01918         if (feasibleSlack>0.0&&newSlack>0.0) {
01919           double smallGap2=smallGap;
01920           if (fabs(0.1*newPrimal)>smallGap) {
01921             smallGap2=0.1*fabs(newPrimal);
01922           } 
01923           if (!fakeUpper(iColumn)) {
01924             double larger;
01925             if (newSlack>feasibleSlack) {
01926               larger=newSlack;
01927             } else {
01928               larger=feasibleSlack;
01929             } 
01930             if (fabs(feasibleSlack-newSlack)<1.0e-6*larger) {
01931               newSlack=feasibleSlack;
01932             } 
01933             //set FAKE here
01934             if (newSlack>wValue2*largestRatio&&newSlack>smallGap2) {
01935               setFakeUpper(iColumn);
01936               newSlack=wValue2*largestRatio;
01937               if (newSlack<smallGap2) {
01938                 newSlack=smallGap2;
01939               } 
01940               numberDecreased++;
01941             } 
01942           } else {
01943             newSlack=wValue2*largestRatio;
01944             if (newSlack<smallGap2) {
01945               newSlack=smallGap2;
01946             } 
01947             if (newSlack>largeGap) {
01948               //increase up to smaller of w.. and largeGap
01949               newSlack=largeGap;
01950             } 
01951             if (newSlack>widenGap*oldSlack) {
01952               numberIncreased++;
01953               //cout<<"wider "<<-j<<" "<<newSlack<<" "<<oldSlack<<" "<<feasibleSlack<<endl;
01954             } 
01955             if (newSlack>feasibleSlack) {
01956               newSlack=feasibleSlack;
01957               clearFakeUpper(iColumn);
01958             } 
01959           } 
01960         } 
01961         if (wVec[iColumn]>dualTolerance) {
01962           dualObjectiveThis-=upper[iColumn]*wVec[iColumn];
01963         } 
01964         upperSlack[iColumn]=newSlack;
01965         if (newSlack<smallerSlack) {
01966           smallerSlack=newSlack;
01967         } 
01968         if (!fakeUpper(iColumn)) {
01969           double infeasibility = fabs(newPrimal+upperSlack[iColumn]-upper[iColumn]);
01970           if (infeasibility>maximumBoundInfeasibility) {
01971             maximumBoundInfeasibility=infeasibility;
01972           } 
01973         } 
01974         if (upperSlack[iColumn]<=kill&&fabs(newPrimal-upper[iColumn])<=kill) {
01975           //may be better to leave at value?
01976           newPrimal=upper[iColumn];
01977           upperSlack[iColumn]=0.0;
01978           thisKilled=true;
01979         } else {
01980           t+=upperSlack[iColumn];
01981         } 
01982       } 
01983       if ((!lowerBound(iColumn))&&
01984                (!upperBound(iColumn))) {
01985         double gap;
01986         if (fabs(newPrimal)>eFree) {
01987           gap=fabs(newPrimal);
01988         } else {
01989           gap=eFree;
01990         } 
01991         zVec[iColumn]=gap*freeMultiplier;
01992         wVec[iColumn]=zVec[iColumn];
01993         //fake to give correct result
01994         s=1.0;
01995         t=1.0;
01996         wValue=0.0;
01997         zValue=2.0*freeMultiplier+qDiagonal;
01998       } 
01999       primal[iColumn]=newPrimal;
02000       if (fabs(newPrimal)>norm) {
02001         norm=fabs(newPrimal);
02002       } 
02003       if (!thisKilled) {
02004         primalObjectiveValue+=primalObjectiveThis;
02005         dualObjectiveValue+=dualObjectiveThis;
02006         if (s>largeGap) {
02007           s=largeGap;
02008         } 
02009         if (t>largeGap) {
02010           t=largeGap;
02011         } 
02012         double divisor = s*wValue+t*zValue;
02013         double diagonalValue=(t*s)/divisor;
02014         diagonal[iColumn]=diagonalValue;
02015         //FUDGE
02016         if (diagonalValue>diagonalLimit) {
02017           //cout<<"large diagonal "<<diagonalValue<<endl;
02018           diagonal[iColumn]=diagonalLimit;
02019         } 
02020         if (diagonalValue<1.0e-10) {
02021           //cout<<"small diagonal "<<diagonalValue<<endl;
02022         } 
02023         if (diagonalValue>largestDiagonal) {
02024           largestDiagonal=diagonalValue;
02025         } 
02026         if (diagonalValue<smallestDiagonal) {
02027           smallestDiagonal=diagonalValue;
02028         } 
02029         double dualInfeasibility=reducedCost-zVec[iColumn]+wVec[iColumn];
02030         if (fabs(dualInfeasibility)>dualTolerance) {
02031           if ((!lowerBound(iColumn))&&
02032                  (!upperBound(iColumn))) {
02033             dualFake+=newPrimal*dualInfeasibility;
02034           } else {
02035             dualFake+=newPrimal*dualInfeasibility;
02036           } 
02037         } 
02038         dualInfeasibility=fabs(dualInfeasibility);
02039         if (dualInfeasibility>maximumDualError) {
02040           maximumDualError=dualInfeasibility;
02041         } 
02042         work1[iColumn]=0.0;
02043       } else {
02044         numberKilled++;
02045         diagonal[iColumn]=0.0;
02046         zVec[iColumn]=0.0;
02047         wVec[iColumn]=0.0;
02048         setFlagged(iColumn);
02049         setFixedOrFree(iColumn);
02050         work1[iColumn]=newPrimal;
02051         offsetObjective+=newPrimal*cost[iColumn];
02052       } 
02053     } else {
02054       work1[iColumn]=primal[iColumn];
02055       diagonal[iColumn]=0.0;
02056       offsetObjective+=primal[iColumn]*cost[iColumn];
02057       if (upper[iColumn]-lower[iColumn]>1.0e-5) {
02058         if (primal[iColumn]<lower[iColumn]+1.0e-8&&dual[iColumn]<-1.0e-8) {
02059           if (-dual[iColumn]>maximumDJInfeasibility) {
02060             maximumDJInfeasibility=-dual[iColumn];
02061           } 
02062         } 
02063         if (primal[iColumn]>upper[iColumn]-1.0e-8&&dual[iColumn]>1.0e-8) {
02064           if (dual[iColumn]>maximumDJInfeasibility) {
02065             maximumDJInfeasibility=dual[iColumn];
02066           } 
02067         } 
02068       } 
02069     } 
02070   } 
02071   handler_->message(CLP_BARRIER_DIAGONAL,messages_)
02072     <<largestDiagonal<<smallestDiagonal
02073     <<CoinMessageEol;
02074   // update rhs region
02075   multiplyAdd(work1+numberColumns_,numberRows_,-1.0,rhsFixRegion_,1.0);
02076   matrix_->times(1.0,work1,rhsFixRegion_);
02077   primalObjectiveValue=primalObjectiveValue+offsetObjective;
02078   dualObjectiveValue+=offsetObjective+dualFake;
02079   if (numberIncreased||numberDecreased) {
02080     handler_->message(CLP_BARRIER_SLACKS,messages_)
02081       <<numberIncreased<<numberDecreased
02082       <<CoinMessageEol;
02083   } 
02084   if (maximumDJInfeasibility) {
02085     handler_->message(CLP_BARRIER_DUALINF,messages_)
02086       <<maximumDJInfeasibility
02087       <<CoinMessageEol;
02088   } 
02089   // Need to rethink (but it is only for printing)
02090   sumPrimalInfeasibilities_=maximumRhsInfeasibility;
02091   sumDualInfeasibilities_=maximumDualError;
02092   maximumBoundInfeasibility_ = maximumBoundInfeasibility;
02093   solutionNorm_ = norm;
02094   //compute error and fixed RHS
02095   multiplyAdd(solution_+numberColumns_,numberRows_,-1.0,errorRegion_,1.0);
02096   matrix_->times(1.0,solution_,errorRegion_);
02097   maximumDualError_=maximumDualError;
02098   maximumBoundInfeasibility_=maximumBoundInfeasibility;
02099   solutionNorm_=solutionNorm;
02100   //finish off objective computation
02101   primalObjective_=primalObjectiveValue*scaleFactor_;
02102   double dualValue2=innerProduct(dual_,numberRows_,
02103                 rhsFixRegion_);
02104   dualObjectiveValue-=dualValue2;
02105   dualObjective_=dualObjectiveValue*scaleFactor_;
02106   if (numberKilled) {
02107       handler_->message(CLP_BARRIER_KILLED,messages_)
02108         <<numberKilled
02109         <<CoinMessageEol;
02110   } 
02111   double maximumRHSError1=0.0;
02112   double maximumRHSError2=0.0;
02113   double primalOffset=0.0;
02114   char * dropped = cholesky_->rowsDropped();
02115   int iRow;
02116   for (iRow=0;iRow<numberRows_;iRow++) {
02117     double value=errorRegion_[iRow];
02118     if (!dropped[iRow]) {
02119       if (fabs(value)>maximumRHSError1) {
02120         maximumRHSError1=fabs(value);
02121       } 
02122     } else {
02123       if (fabs(value)>maximumRHSError2) {
02124         maximumRHSError2=fabs(value);
02125       } 
02126       primalOffset+=value*dual_[iRow];
02127     } 
02128   } 
02129   primalObjective_-=primalOffset*scaleFactor_;
02130   if (maximumRHSError1>maximumRHSError2) {
02131     maximumRHSError_=maximumRHSError1;
02132   } else {
02133     maximumRHSError_=maximumRHSError1; //note change
02134     if (maximumRHSError2>primalTolerance()) {
02135       handler_->message(CLP_BARRIER_ABS_DROPPED,messages_)
02136         <<maximumRHSError2
02137         <<CoinMessageEol;
02138     } 
02139   } 
02140   objectiveNorm_=maximumAbsElement(dual_,numberRows_);
02141   if (objectiveNorm_<1.0e-12) {
02142     objectiveNorm_=1.0e-12;
02143   } 
02144   if (objectiveNorm_<baseObjectiveNorm_) {
02145     //std::cout<<" base "<<baseObjectiveNorm_<<" "<<objectiveNorm_<<std::endl;
02146     if (objectiveNorm_<baseObjectiveNorm_*1.0e-4) {
02147       objectiveNorm_=baseObjectiveNorm_*1.0e-4;
02148     } 
02149   } 
02150   bool primalFeasible=true;
02151   if (maximumRHSError_>primalTolerance()||
02152       maximumDualError_>dualTolerance/scaleFactor_) {
02153     handler_->message(CLP_BARRIER_ABS_ERROR,messages_)
02154       <<maximumRHSError_<<maximumDualError_
02155       <<CoinMessageEol;
02156   } 
02157   if (rhsNorm_>solutionNorm_) {
02158     solutionNorm_=rhsNorm_;
02159   } 
02160   double scaledRHSError=maximumRHSError_/solutionNorm_;
02161   bool dualFeasible=true;
02162   if (maximumBoundInfeasibility_>primalTolerance()||
02163       scaledRHSError>primalTolerance())
02164     primalFeasible=false;
02165   if (maximumDualError_>objectiveNorm_*dualTolerance) 
02166     dualFeasible=false;
02167   if (!primalFeasible||!dualFeasible) {
02168     handler_->message(CLP_BARRIER_FEASIBLE,messages_)
02169       <<maximumBoundInfeasibility_<<scaledRHSError
02170       <<maximumDualError_/objectiveNorm_
02171       <<CoinMessageEol;
02172   } 
02173   if (!gonePrimalFeasible_) {
02174     gonePrimalFeasible_=primalFeasible;
02175   } else if (!primalFeasible) {
02176     gonePrimalFeasible_=primalFeasible;
02177     if (!numberKilled) {
02178       handler_->message(CLP_BARRIER_GONE_INFEASIBLE,messages_)
02179       <<CoinMessageEol;
02180     } 
02181   } 
02182   if (!goneDualFeasible_) {
02183     goneDualFeasible_=dualFeasible;
02184   } else if (!dualFeasible) {
02185     handler_->message(CLP_BARRIER_GONE_INFEASIBLE,messages_)
02186       <<CoinMessageEol;
02187     goneDualFeasible_=dualFeasible;
02188   } 
02189   //objectiveValue();
02190   if (solutionNorm_>1.0e40) {
02191     std::cout <<"primal off to infinity"<<std::endl;
02192     abort();
02193   } 
02194   if (objectiveNorm_>1.0e40) {
02195     std::cout <<"dual off to infinity"<<std::endl;
02196     abort();
02197   } 
02198   handler_->message(CLP_BARRIER_STEP,messages_)
02199     <<actualPrimalStep_
02200     <<actualDualStep_
02201     <<mu_
02202     <<CoinMessageEol;
02203   numberIterations_++;
02204   return numberKilled;
02205 }
02206 //  Save info on products of affine deltaT*deltaW and deltaS*deltaZ
02207 double 
02208 ClpPredictorCorrector::affineProduct()
02209 {
02210   double * work1 = deltaZ_;
02211   double * work2 = deltaW_;
02212   double * work3 = deltaS_;
02213   double * work4 = deltaT_;
02214   double * lower = lower_;
02215   double * upper = upper_;
02216   double * lowerSlack = lowerSlack_;
02217   double * upperSlack = upperSlack_;
02218   double * primal = solution_;
02219   //direction vector in weights
02220   double * weights = weights_;
02221   double product = 0.0;
02222   //IF zVec starts as 0 then work1 always zero
02223   //(remember if free then zVec not 0)
02224   //I think free can be done with careful use of boundSlacks to zero
02225   //out all we want
02226   for (int iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
02227     double w3=work1[iColumn]*weights[iColumn];
02228     double w4=-work2[iColumn]*weights[iColumn];
02229     if (lowerBound(iColumn)) {
02230       if (!fakeLower(iColumn)) {
02231         w3+=work1[iColumn]*(primal[iColumn]-lowerSlack[iColumn]-lower[iColumn]);
02232       } 
02233       product+=w3;
02234     } 
02235     if (upperBound(iColumn)) {
02236       if (!fakeUpper(iColumn)) {
02237         w4+=work2[iColumn]*(-primal[iColumn]-upperSlack[iColumn]+upper[iColumn]);
02238       } 
02239       product+=w4;
02240     } 
02241     work3[iColumn]=w3;
02242     work4[iColumn]=w4;
02243   } 
02244   return product;
02245 }