source: trunk/Clp/src/ClpPredictorCorrector.cpp @ 1665

Last change on this file since 1665 was 1665, checked in by lou, 10 years ago

Add EPL license notice in src.

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