source: stable/1.15/Clp/src/ClpPredictorCorrector.cpp @ 2018

Last change on this file since 2018 was 1959, checked in by stefan, 6 years ago

merge r1945, r1947, and r1950..r1958 from trunk

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