source: stable/1.17/Clp/src/ClpPredictorCorrector.cpp @ 2439

Last change on this file since 2439 was 2385, checked in by unxusr, 11 months ago

formatting

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