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

Last change on this file since 1926 was 1926, checked in by forrest, 7 years ago

try to fix infeasibility ray,
changes as in stable (for presolve),
stuff for Cbc parameters

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