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

Last change on this file since 1912 was 1726, checked in by stefan, 9 years ago

fix compiler warnings, including one that pointed to a bug

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