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

Last change on this file since 1949 was 1931, checked in by stefan, 7 years ago

sync with trunk rev 1930

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