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

Last change on this file since 1577 was 1525, checked in by mjs, 10 years ago

Formatted .cpp, .hpp, .c, .h files with "astyle -A4 -p". This matches the formatting used in the grand CBC reorganization.

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