# source:trunk/Clp/src/ClpSimplexPrimal.cpp@2447

Last change on this file since 2447 was 2447, checked in by forrest, 7 months ago

• Property svn:eol-style set to `native`
• Property svn:keywords set to `Id`
File size: 146.9 KB
Line
1/* \$Id: ClpSimplexPrimal.cpp 2447 2019-04-07 09:56:25Z forrest \$ */
5
6/* Notes on implementation of primal simplex algorithm.
7
8   When primal feasible(A):
9
10   If dual feasible, we are optimal.  Otherwise choose an infeasible
11   basic variable to enter basis from a bound (B).  We now need to find an
12   outgoing variable which will leave problem primal feasible so we get
13   the column of the tableau corresponding to the incoming variable
14   (with the correct sign depending if variable will go up or down).
15
16   We now perform a ratio test to determine which outgoing variable will
17   preserve primal feasibility (C).  If no variable found then problem
18   is unbounded (in primal sense).  If there is a variable, we then
19   perform pivot and repeat.  Trivial?
20
21   -------------------------------------------
22
23   A) How do we get primal feasible?  All variables have fake costs
24   outside their feasible region so it is trivial to declare problem
25   feasible.  OSL did not have a phase 1/phase 2 approach but
26   instead effectively put an extra cost on infeasible basic variables
27   I am taking the same approach here, although it is generalized
28   to allow for non-linear costs and dual information.
29
30   In OSL, this weight was changed heuristically, here at present
31   it is only increased if problem looks finished.  If problem is
32   feasible I check for unboundedness.  If not unbounded we
33   could play with going into dual.  As long as weights increase
34   any algorithm would be finite.
35
36   B) Which incoming variable to choose is a virtual base class.
37   For difficult problems steepest edge is preferred while for
38   very easy (large) problems we will need partial scan.
39
40   C) Sounds easy, but this is hardest part of algorithm.
41      1) Instead of stopping at first choice, we may be able
42      to allow that variable to go through bound and if objective
43      still improving choose again.  These mini iterations can
44      increase speed by orders of magnitude but we may need to
45      go to more of a bucket choice of variable rather than looking
46      at them one by one (for speed).
47      2) Accuracy.  Basic infeasibilities may be less than
48      tolerance.  Pivoting on these makes objective go backwards.
49      OSL modified cost so a zero move was made, Gill et al
50      modified so a strictly positive move was made.
51      The two problems are that re-factorizations can
52      change rinfeasibilities above and below tolerances and that when
53      finished we need to reset costs and try again.
54      3) Degeneracy.  Gill et al helps but may not be enough.  We
55      may need more.  Also it can improve speed a lot if we perturb
56      the rhs and bounds significantly.
57
58  References:
59     Forrest and Goldfarb, Steepest-edge simplex algorithms for
60       linear programming - Mathematical Programming 1992
61     Forrest and Tomlin, Implementing the simplex method for
62       the Optimization Subroutine Library - IBM Systems Journal 1992
63     Gill, Murray, Saunders, Wright A Practical Anti-Cycling
64       Procedure for Linear and Nonlinear Programming SOL report 1988
65
66
67  TODO:
68
69  a) Better recovery procedures.  At present I never check on forward
71     re-factorization frequency, but this is only on singular
72     factorizations.
73  b) Fast methods for large easy problems (and also the option for
74     the code to automatically choose which method).
75  c) We need to be able to stop in various ways for OSI - this
76     is fairly easy.
77
78 */
79
80#include "CoinPragma.hpp"
81
82#include <math.h>
83//#define FAKE_CILK
84
85#include "CoinHelperFunctions.hpp"
86#include "ClpSimplexPrimal.hpp"
87#include "ClpFactorization.hpp"
88#include "ClpNonLinearCost.hpp"
89#include "CoinPackedMatrix.hpp"
90#include "CoinIndexedVector.hpp"
91#include "ClpPrimalColumnPivot.hpp"
92#include "ClpMessage.hpp"
93#include "ClpEventHandler.hpp"
94#include <cfloat>
95#include <cassert>
96#include <string>
97#include <stdio.h>
98#include <iostream>
99#if ABOCA_LITE
100#undef ALT_UPDATE_WEIGHTS
101#define ALT_UPDATE_WEIGHTS 2
102#endif
103#if ALT_UPDATE_WEIGHTS == 1
104CoinIndexedVector *altVector[3] = { NULL, NULL, NULL };
105#endif
106#ifdef CLP_USER_DRIVEN1
107/* Returns true if variable sequenceOut can leave basis when
108   model->sequenceIn() enters.
109   This function may be entered several times for each sequenceOut.
110   The first time realAlpha will be positive if going to lower bound
111   and negative if going to upper bound (scaled bounds in lower,upper) - then will be zero.
112   currentValue is distance to bound.
113   currentTheta is current theta.
114   alpha is fabs(pivot element).
115   Variable will change theta if currentValue - currentTheta*alpha < 0.0
116*/
117bool userChoiceValid1(const ClpSimplex *model,
118  int sequenceOut,
119  double currentValue,
120  double currentTheta,
121  double alpha,
122  double realAlpha);
123/* This returns true if chosen in/out pair valid.
124   The main thing to check would be variable flipping bounds may be
125   OK.  This would be signaled by reasonable theta_ and valueOut_.
126   If you return false sequenceIn_ will be flagged as ineligible.
127*/
128bool userChoiceValid2(const ClpSimplex *model);
129/* If a good pivot then you may wish to unflag some variables.
130 */
131void userChoiceWasGood(ClpSimplex *model);
132#endif
133// primal
134int ClpSimplexPrimal::primal(int ifValuesPass, int startFinishOptions)
135{
136
137  /*
138         Method
139
140        It tries to be a single phase approach with a weight of 1.0 being
141        given to getting optimal and a weight of infeasibilityCost_ being
142        given to getting primal feasible.  In this version I have tried to
143        be clever in a stupid way.  The idea of fake bounds in dual
144        seems to work so the primal analogue would be that of getting
145        bounds on reduced costs (by a presolve approach) and using
146        these for being above or below feasible region.  I decided to waste
147        memory and keep these explicitly.  This allows for non-linear
148        costs!
149
150        The code is designed to take advantage of sparsity so arrays are
151        seldom zeroed out from scratch or gone over in their entirety.
152        The only exception is a full scan to find incoming variable for
153        Dantzig row choice.  For steepest edge we keep an updated list
154        of dual infeasibilities (actually squares).
155        On easy problems we don't need full scan - just
156        pick first reasonable.
157
158        One problem is how to tackle degeneracy and accuracy.  At present
159        I am using the modification of costs which I put in OSL and which was
160        extended by Gill et al.  I am still not sure of the exact details.
161
162        The flow of primal is three while loops as follows:
163
164        while (not finished) {
165
166          while (not clean solution) {
167
168             Factorize and/or clean up solution by changing bounds so
169       primal feasible.  If looks finished check fake primal bounds.
170       Repeat until status is iterating (-1) or finished (0,1,2)
171
172          }
173
174          while (status==-1) {
175
176            Iterate until no pivot in or out or time to re-factorize.
177
178            Flow is:
179
180            choose pivot column (incoming variable).  if none then
181      we are primal feasible so looks as if done but we need to
182      break and check bounds etc.
183
184      Get pivot column in tableau
185
186            Choose outgoing row.  If we don't find one then we look
187      primal unbounded so break and check bounds etc.  (Also the
188      pivot tolerance is larger after any iterations so that may be
189      reason)
190
191            If we do find outgoing row, we may have to adjust costs to
192      keep going forwards (anti-degeneracy).  Check pivot will be stable
193      and if unstable throw away iteration and break to re-factorize.
194      If minor error re-factorize after iteration.
195
196      Update everything (this may involve changing bounds on
197      variables to stay primal feasible.
198
199          }
200
201        }
202
203        At present we never check we are going forwards.  I overdid that in
204        OSL so will try and make a last resort.
205
206        Needs partial scan pivot in option.
207
208        May need other anti-degeneracy measures, especially if we try and use
209        loose tolerances as a way to solve in fewer iterations.
210
211        I like idea of dynamic scaling.  This gives opportunity to decouple
212        different implications of scaling for accuracy, iteration count and
213        feasibility tolerance.
214
215     */
216
217  algorithm_ = +1;
218  moreSpecialOptions_ &= ~16; // clear check replaceColumn accuracy
219
220  // save data
221  ClpDataSave data = saveData();
222  if (problemStatus_ == 10 && sumPrimalInfeasibilities_ == -123456789.0) {
223    // large infeasibility cost wanted
224    infeasibilityCost_ = CoinMax(infeasibilityCost_, 1.0e13);
225  }
226  matrix_->refresh(this); // make sure matrix okay
227
228  // Save so can see if doing after dual
229  int initialStatus = problemStatus_;
230  int initialIterations = numberIterations_;
231  int initialNegDjs = -1;
232  // initialize - maybe values pass and algorithm_ is +1
233#if 0
234     // if so - put in any superbasic costed slacks
235     if (ifValuesPass && specialOptions_ < 0x01000000) {
236          // Get column copy
237          const CoinPackedMatrix * columnCopy = matrix();
238          const int * row = columnCopy->getIndices();
239          const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
240          const int * columnLength = columnCopy->getVectorLengths();
241          //const double * element = columnCopy->getElements();
242          int n = 0;
243          for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
244               if (columnLength[iColumn] == 1) {
245                    Status status = getColumnStatus(iColumn);
246                    if (status != basic && status != isFree) {
247                         double value = columnActivity_[iColumn];
248                         if (fabs(value - columnLower_[iColumn]) > primalTolerance_ &&
249                                   fabs(value - columnUpper_[iColumn]) > primalTolerance_) {
250                              int iRow = row[columnStart[iColumn]];
251                              if (getRowStatus(iRow) == basic) {
252                                   setRowStatus(iRow, superBasic);
253                                   setColumnStatus(iColumn, basic);
254                                   n++;
255                              }
256                         }
257                    }
258               }
259          }
260          printf("%d costed slacks put in basis\n", n);
261     }
262#endif
263  // Start can skip some things in transposeTimes
264  specialOptions_ |= 131072;
265  if (!startup(ifValuesPass, startFinishOptions)) {
266    // See if better to use all slack basis
267    if (nonLinearCost_->sumInfeasibilities() > 1.0e15) {
268      // If not all slack then make it
269      int numberRowBasic = 0;
270      for (int i = 0; i < numberRows_; i++) {
271        if (getRowStatus(i) == basic)
272          numberRowBasic++;
273      }
274      if (numberRowBasic < numberRows_) {
275        allSlackBasis(true);
276        int lastCleaned = -10000;
277        statusOfProblemInPrimal(lastCleaned, 1, &progress_, true, ifValuesPass, NULL);
278      }
279    }
280
281    // Set average theta
282    nonLinearCost_->setAverageTheta(1.0e3);
283    int lastCleaned = 0; // last time objective or bounds cleaned up
284
285    // Say no pivot has occurred (for steepest edge and updates)
286    pivotRow_ = -2;
287
288    // This says whether to restore things etc
289    int factorType = 0;
290    if (problemStatus_ < 0 && perturbation_ < 100 && !ifValuesPass) {
291      perturb(0);
292      // Can't get here if values pass
293      assert(!ifValuesPass);
294      gutsOfSolution(NULL, NULL);
295      if (handler_->logLevel() > 2) {
296        handler_->message(CLP_SIMPLEX_STATUS, messages_)
297          << numberIterations_ << objectiveValue();
298        handler_->printing(sumPrimalInfeasibilities_ > 0.0)
299          << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
300        handler_->printing(sumDualInfeasibilities_ > 0.0)
301          << sumDualInfeasibilities_ << numberDualInfeasibilities_;
302        handler_->printing(numberDualInfeasibilitiesWithoutFree_
303          < numberDualInfeasibilities_)
304          << numberDualInfeasibilitiesWithoutFree_;
305        handler_->message() << CoinMessageEol;
306      }
307    }
308    ClpSimplex *saveModel = NULL;
309    int stopSprint = -1;
310    int sprintPass = 0;
311    int reasonableSprintIteration = 0;
312    int lastSprintIteration = 0;
313    double lastObjectiveValue = COIN_DBL_MAX;
314    // Start check for cycles
315    progress_.fillFromModel(this);
316    progress_.startCheck();
317    /*
318            Status of problem:
319            0 - optimal
320            1 - infeasible
321            2 - unbounded
322            -1 - iterating
323            -2 - factorization wanted
324            -3 - redo checking without factorization
325            -4 - looks infeasible
326            -5 - looks unbounded
327          */
328    while (problemStatus_ < 0) {
329      int iRow, iColumn;
330      // clear
331      for (iRow = 0; iRow < 4; iRow++) {
332        rowArray_[iRow]->clear();
333      }
334
335      for (iColumn = 0; iColumn < SHORT_REGION; iColumn++) {
336        columnArray_[iColumn]->clear();
337      }
338
339      // give matrix (and model costs and bounds a chance to be
340      // refreshed (normally null)
341      matrix_->refresh(this);
342      // If getting nowhere - why not give it a kick
343#if 1
344      if (perturbation_ < 101 && numberIterations_ > 2 * (numberRows_ + numberColumns_) && (specialOptions_ & 4) == 0
345        && initialStatus != 10) {
346        perturb(1);
347        matrix_->rhsOffset(this, true, false);
348      }
349#endif
350      // If we have done no iterations - special
351      if (lastGoodIteration_ == numberIterations_ && factorType)
352        factorType = 3;
353      if (saveModel) {
354        // Doing sprint
355        if (sequenceIn_ < 0 || numberIterations_ >= stopSprint) {
356          problemStatus_ = -1;
357          originalModel(saveModel);
358          saveModel = NULL;
359          if (sequenceIn_ < 0 && numberIterations_ < reasonableSprintIteration && sprintPass > 100)
360            primalColumnPivot_->switchOffSprint();
361          //lastSprintIteration=numberIterations_;
362          COIN_DETAIL_PRINT(printf("End small model\n"));
363        }
364      }
365
366      // may factorize, checks if problem finished
367      statusOfProblemInPrimal(lastCleaned, factorType, &progress_, true, ifValuesPass, saveModel);
368      if (initialStatus == 10) {
369        // cleanup phase
370        if (initialIterations != numberIterations_) {
371          if (numberDualInfeasibilities_ > 10000 && numberDualInfeasibilities_ > 10 * initialNegDjs) {
372            // getting worse - try perturbing
373            if (perturbation_ < 101 && (specialOptions_ & 4) == 0) {
374              perturb(1);
375              matrix_->rhsOffset(this, true, false);
376              statusOfProblemInPrimal(lastCleaned, factorType, &progress_, true, ifValuesPass, saveModel);
377            }
378          }
379        } else {
380          // save number of negative djs
381          if (!numberPrimalInfeasibilities_)
382            initialNegDjs = numberDualInfeasibilities_;
383          // make sure weight won't be changed
384          if (infeasibilityCost_ == 1.0e10)
385            infeasibilityCost_ = 1.000001e10;
386        }
387      }
388      // See if sprint says redo because of problems
389      if (numberDualInfeasibilities_ == -776) {
390        // Need new set of variables
391        problemStatus_ = -1;
392        originalModel(saveModel);
393        saveModel = NULL;
394        //lastSprintIteration=numberIterations_;
395        COIN_DETAIL_PRINT(printf("End small model after\n"));
396        statusOfProblemInPrimal(lastCleaned, factorType, &progress_, true, ifValuesPass, saveModel);
397      }
398      int numberSprintIterations = 0;
399      int numberSprintColumns = primalColumnPivot_->numberSprintColumns(numberSprintIterations);
400      if (problemStatus_ == 777) {
401        // problems so do one pass with normal
402        problemStatus_ = -1;
403        originalModel(saveModel);
404        saveModel = NULL;
405        // Skip factorization
406        //statusOfProblemInPrimal(lastCleaned,factorType,&progress_,false,saveModel);
407        statusOfProblemInPrimal(lastCleaned, factorType, &progress_, true, ifValuesPass, saveModel);
408      } else if (problemStatus_ < 0 && !saveModel && numberSprintColumns && firstFree_ < 0) {
409        int numberSort = 0;
410        int numberFixed = 0;
411        int numberBasic = 0;
412        reasonableSprintIteration = numberIterations_ + 100;
413        int *whichColumns = new int[numberColumns_];
414        double *weight = new double[numberColumns_];
415        int numberNegative = 0;
416        double sumNegative = 0.0;
417        // now massage weight so all basic in plus good djs
418        for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
419          double dj = dj_[iColumn];
420          switch (getColumnStatus(iColumn)) {
421
422          case basic:
423            dj = -1.0e50;
424            numberBasic++;
425            break;
426          case atUpperBound:
427            dj = -dj;
428            break;
429          case isFixed:
430            dj = 1.0e50;
431            numberFixed++;
432            break;
433          case atLowerBound:
434            dj = dj;
435            break;
436          case isFree:
437            dj = -100.0 * fabs(dj);
438            break;
439          case superBasic:
440            dj = -100.0 * fabs(dj);
441            break;
442          }
443          if (dj < -dualTolerance_ && dj > -1.0e50) {
444            numberNegative++;
445            sumNegative -= dj;
446          }
447          weight[iColumn] = dj;
448          whichColumns[iColumn] = iColumn;
449        }
450        handler_->message(CLP_SPRINT, messages_)
451          << sprintPass << numberIterations_ - lastSprintIteration << objectiveValue() << sumNegative
452          << numberNegative
453          << CoinMessageEol;
454        sprintPass++;
455        lastSprintIteration = numberIterations_;
456        if (objectiveValue() * optimizationDirection_ > lastObjectiveValue - 1.0e-7 && sprintPass > 5) {
457          // switch off
458          COIN_DETAIL_PRINT(printf("Switching off sprint\n"));
459          primalColumnPivot_->switchOffSprint();
460        } else {
461          lastObjectiveValue = objectiveValue() * optimizationDirection_;
462          // sort
463          CoinSort_2(weight, weight + numberColumns_, whichColumns);
464          numberSort = CoinMin(numberColumns_ - numberFixed, numberBasic + numberSprintColumns);
465          // Sort to make consistent ?
466          std::sort(whichColumns, whichColumns + numberSort);
467          saveModel = new ClpSimplex(this, numberSort, whichColumns);
468          delete[] whichColumns;
469          delete[] weight;
470          // Skip factorization
471          //statusOfProblemInPrimal(lastCleaned,factorType,&progress_,false,saveModel);
472          //statusOfProblemInPrimal(lastCleaned,factorType,&progress_,true,saveModel);
473          stopSprint = numberIterations_ + numberSprintIterations;
474          COIN_DETAIL_PRINT(printf("Sprint with %d columns for %d iterations\n",
475            numberSprintColumns, numberSprintIterations));
476        }
477      }
478
479      // Say good factorization
480      factorType = 1;
481
482      // Say no pivot has occurred (for steepest edge and updates)
483      pivotRow_ = -2;
484      // exit if feasible and done enough iterations
485      if ((moreSpecialOptions_ & 1048576) != 0) {
486        int maxIterations = maximumIterations() - 1000000;
487        if (maxIterations > 0 && maxIterations < 200000) {
488          if (!nonLinearCost_->numberInfeasibilities() && numberIterations_ >= maxIterations) {
489            problemStatus_ = 3;
490            secondaryStatus_ = 10;
491          }
492        }
493      }
494      // exit if victory declared
495      if (problemStatus_ >= 0) {
496#ifdef CLP_USER_DRIVEN
497        int status = eventHandler_->event(ClpEventHandler::endInPrimal);
498        if (status >= 0 && status < 10) {
499          // carry on
500          problemStatus_ = -1;
501          if (status == 0)
502            continue; // re-factorize
503        } else if (status >= 10) {
504          problemStatus_ = status - 10;
505          break;
506        } else {
507          break;
508        }
509#else
510        break;
511#endif
512      }
513
514      // test for maximum iterations
515      if (hitMaximumIterations() || (ifValuesPass == 2 && firstFree_ < 0)) {
516        problemStatus_ = 3;
517        break;
518      } else if ((moreSpecialOptions_ & 524288) != 0 && !nonLinearCost_->numberInfeasibilities() && fabs(dblParam_[ClpDualObjectiveLimit]) > 1.0e30) {
519        problemStatus_ = 3;
520        secondaryStatus_ = 10;
521        break;
522      }
523
524      if (firstFree_ < 0) {
525        if (ifValuesPass) {
526          // end of values pass
527          ifValuesPass = 0;
528          int status = eventHandler_->event(ClpEventHandler::endOfValuesPass);
529          if (status >= 0) {
530            problemStatus_ = 5;
531            secondaryStatus_ = ClpEventHandler::endOfValuesPass;
532            break;
533          }
534          //#define FEB_TRY
535#if 1 //def FEB_TRY
536          if (perturbation_ < 100)
537            perturb(0);
538#endif
539        }
540      }
541      // Check event
542      {
543        int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
544        if (status >= 0) {
545          // if >=100 - then special e.g. unperturb
546          if (status != 101) {
547            problemStatus_ = 5;
548            secondaryStatus_ = ClpEventHandler::endOfFactorization;
549            break;
550          } else {
551            unPerturb();
552            continue;
553          }
554        }
555      }
556      // Iterate
557      whileIterating(ifValuesPass ? 1 : 0);
558      if (sequenceIn_ < 0 && ifValuesPass == 2)
559        problemStatus_ = 3; // user wants to exit
560    }
561  }
562  // if infeasible get real values
563  //printf("XXXXY final cost %g\n",infeasibilityCost_);
564  progress_.initialWeight_ = 0.0;
565  if (problemStatus_ == 1 && secondaryStatus_ != 6) {
566    double saveWeight = infeasibilityCost_;
567#ifndef WANT_INFEASIBLE_DUALS
568    infeasibilityCost_ = 0.0;
569    createRim(1 + 4);
570#else
571    infeasibilityCost_ = 1.0;
572    createRim(1);
573    memset(cost_, 0, (numberRows_ + numberColumns_) * sizeof(double));
574#endif
575    delete nonLinearCost_;
576    nonLinearCost_ = new ClpNonLinearCost(this);
577    nonLinearCost_->checkInfeasibilities(0.0);
578    sumPrimalInfeasibilities_ = nonLinearCost_->sumInfeasibilities();
579    numberPrimalInfeasibilities_ = nonLinearCost_->numberInfeasibilities();
580    // and get good feasible duals
581    computeDuals(NULL);
582    infeasibilityCost_ = saveWeight;
583  }
584  // Stop can skip some things in transposeTimes
585  specialOptions_ &= ~131072;
586  // clean up
587  unflag();
588  finish(startFinishOptions);
589  restoreData(data);
590  return problemStatus_;
591}
592/*
593  Reasons to come out:
594  -1 iterations etc
595  -2 inaccuracy
596  -3 slight inaccuracy (and done iterations)
597  -4 end of values pass and done iterations
598  +0 looks optimal (might be infeasible - but we will investigate)
599  +2 looks unbounded
600  +3 max iterations
601*/
602int ClpSimplexPrimal::whileIterating(int valuesOption)
603{
604  // Say if values pass
605  int ifValuesPass = (firstFree_ >= 0) ? 1 : 0;
606  int returnCode = -1;
607  int superBasicType = 1;
608  if (valuesOption > 1)
609    superBasicType = 3;
610  // delete any rays
611  delete[] ray_;
612  ray_ = NULL;
613  // status stays at -1 while iterating, >=0 finished, -2 to invert
614  // status -3 to go to top without an invert
615  while (problemStatus_ == -1) {
616    //#define CLP_DEBUG 1
617#ifdef CLP_DEBUG
618    {
619      int i;
620      // not [1] as has information
621      for (i = 0; i < 4; i++) {
622        if (i != 1)
623          rowArray_[i]->checkClear();
624      }
625      for (i = 0; i < 2; i++) {
626        columnArray_[i]->checkClear();
627      }
628    }
629#endif
630#if 0
631          {
632               int iPivot;
633               double * array = rowArray_[3]->denseVector();
634               int * index = rowArray_[3]->getIndices();
635               int i;
636               for (iPivot = 0; iPivot < numberRows_; iPivot++) {
637                    int iSequence = pivotVariable_[iPivot];
638                    unpackPacked(rowArray_[3], iSequence);
639                    factorization_->updateColumn(rowArray_[2], rowArray_[3]);
640                    int number = rowArray_[3]->getNumElements();
641                    for (i = 0; i < number; i++) {
642                         int iRow = index[i];
643                         if (iRow == iPivot)
644                              assert (fabs(array[i] - 1.0) < 1.0e-4);
645                         else
646                              assert (fabs(array[i]) < 1.0e-4);
647                    }
648                    rowArray_[3]->clear();
649               }
650          }
651#endif
652#if 0
653          nonLinearCost_->checkInfeasibilities(primalTolerance_);
654          printf("suminf %g number %d\n", nonLinearCost_->sumInfeasibilities(),
655                 nonLinearCost_->numberInfeasibilities());
656#endif
657#if CLP_DEBUG > 2
658    // very expensive
659    if (numberIterations_ > 0 && numberIterations_ < 100 && !ifValuesPass) {
660      handler_->setLogLevel(63);
661      double saveValue = objectiveValue_;
662      double *saveRow1 = new double[numberRows_];
663      double *saveRow2 = new double[numberRows_];
664      CoinMemcpyN(rowReducedCost_, numberRows_, saveRow1);
665      CoinMemcpyN(rowActivityWork_, numberRows_, saveRow2);
666      double *saveColumn1 = new double[numberColumns_];
667      double *saveColumn2 = new double[numberColumns_];
668      CoinMemcpyN(reducedCostWork_, numberColumns_, saveColumn1);
669      CoinMemcpyN(columnActivityWork_, numberColumns_, saveColumn2);
670      gutsOfSolution(NULL, NULL, false);
671      printf("xxx %d old obj %g, recomputed %g, sum primal inf %g\n",
672        numberIterations_,
673        saveValue, objectiveValue_, sumPrimalInfeasibilities_);
674      CoinMemcpyN(saveRow1, numberRows_, rowReducedCost_);
675      CoinMemcpyN(saveRow2, numberRows_, rowActivityWork_);
676      CoinMemcpyN(saveColumn1, numberColumns_, reducedCostWork_);
677      CoinMemcpyN(saveColumn2, numberColumns_, columnActivityWork_);
678      delete[] saveRow1;
679      delete[] saveRow2;
680      delete[] saveColumn1;
681      delete[] saveColumn2;
682      objectiveValue_ = saveValue;
683    }
684#endif
685    if (!ifValuesPass) {
686      // choose column to come in
687      // can use pivotRow_ to update weights
688      // pass in list of cost changes so can do row updates (rowArray_[1])
689      // NOTE rowArray_[0] is used by computeDuals which is a
690      // slow way of getting duals but might be used
691#ifdef LONG_REGION_2
692      primalColumn(rowArray_[1], rowArray_[0], rowArray_[3],
693        columnArray_[0], rowArray_[2]);
694#else
695      primalColumn(rowArray_[1], rowArray_[2], rowArray_[3],
696        columnArray_[0], columnArray_[1]);
697#endif
698    } else {
699      // in values pass
700      int sequenceIn = nextSuperBasic(superBasicType, columnArray_[0]);
701      if (valuesOption > 1)
702        superBasicType = 2;
703      if (sequenceIn < 0) {
704        // end of values pass - initialize weights etc
705        handler_->message(CLP_END_VALUES_PASS, messages_)
706          << numberIterations_;
707        primalColumnPivot_->saveWeights(this, 5);
708        problemStatus_ = -2; // factorize now
709        pivotRow_ = -1; // say no weights update
710        returnCode = -4;
711        // Clean up
712        int i;
713        for (i = 0; i < numberRows_ + numberColumns_; i++) {
714          if (getColumnStatus(i) == atLowerBound || getColumnStatus(i) == isFixed)
715            solution_[i] = lower_[i];
716          else if (getColumnStatus(i) == atUpperBound)
717            solution_[i] = upper_[i];
718        }
719        break;
720      } else {
721        // normal
722        sequenceIn_ = sequenceIn;
723        valueIn_ = solution_[sequenceIn_];
724        lowerIn_ = lower_[sequenceIn_];
725        upperIn_ = upper_[sequenceIn_];
726        dualIn_ = dj_[sequenceIn_];
727      }
728    }
729    pivotRow_ = -1;
730    sequenceOut_ = -1;
731    rowArray_[1]->clear();
732    if (sequenceIn_ >= 0) {
733      // we found a pivot column
734      assert(!flagged(sequenceIn_));
735#ifdef CLP_DEBUG
736      if ((handler_->logLevel() & 32)) {
737        char x = isColumn(sequenceIn_) ? 'C' : 'R';
738        std::cout << "pivot column " << x << sequenceWithin(sequenceIn_) << std::endl;
739      }
740#endif
741#ifdef CLP_DEBUG
742      {
743        int checkSequence = -2077;
744        if (checkSequence >= 0 && checkSequence < numberRows_ + numberColumns_ && !ifValuesPass) {
745          rowArray_[2]->checkClear();
746          rowArray_[3]->checkClear();
747          double *array = rowArray_[3]->denseVector();
748          int *index = rowArray_[3]->getIndices();
749          unpackPacked(rowArray_[3], checkSequence);
750          factorization_->updateColumnForDebug(rowArray_[2], rowArray_[3]);
751          int number = rowArray_[3]->getNumElements();
752          double dualIn = cost_[checkSequence];
753          int i;
754          for (i = 0; i < number; i++) {
755            int iRow = index[i];
756            int iPivot = pivotVariable_[iRow];
757            double alpha = array[i];
758            dualIn -= alpha * cost_[iPivot];
759          }
760          printf("old dj for %d was %g, recomputed %g\n", checkSequence,
761            dj_[checkSequence], dualIn);
762          rowArray_[3]->clear();
763          if (numberIterations_ > 2000)
764            exit(1);
765        }
766      }
767#endif
768      // do second half of iteration
769      returnCode = pivotResult(ifValuesPass);
770      if (returnCode < -1 && returnCode > -5) {
771        problemStatus_ = -2; //
772      } else if (returnCode == -5) {
773        if ((moreSpecialOptions_ & 16) == 0 && factorization_->pivots()) {
774          moreSpecialOptions_ |= 16;
775          problemStatus_ = -2;
776        }
777        // otherwise something flagged - continue;
778      } else if (returnCode == 2) {
779        problemStatus_ = -5; // looks unbounded
780      } else if (returnCode == 4) {
781        problemStatus_ = -2; // looks unbounded but has iterated
782      } else if (returnCode != -1) {
783        assert(returnCode == 3);
784        if (problemStatus_ != 5)
785          problemStatus_ = 3;
786      }
787    } else {
788      // no pivot column
789#ifdef CLP_DEBUG
790      if (handler_->logLevel() & 32)
791        printf("** no column pivot\n");
792#endif
793      if (nonLinearCost_->numberInfeasibilities())
794        problemStatus_ = -4; // might be infeasible
795      // Force to re-factorize early next time
796      int numberPivots = factorization_->pivots();
797      returnCode = 0;
798#ifdef CLP_USER_DRIVEN
799      // If large number of pivots trap later?
800      if (problemStatus_ == -1 && numberPivots < 1000000) {
801        int status = eventHandler_->event(ClpEventHandler::noCandidateInPrimal);
802        if (status >= 0 && status < 10) {
803          // carry on
804          problemStatus_ = -1;
805          if (status == 0)
806            break;
807        } else if (status >= 10) {
808          problemStatus_ = status - 10;
809          break;
810        } else {
811          forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
812          break;
813        }
814      }
815#else
816      forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
817      break;
818#endif
819    }
820  }
821  if (valuesOption > 1)
822    columnArray_[0]->setNumElements(0);
823  return returnCode;
824}
825/* Checks if finished.  Updates status */
826void ClpSimplexPrimal::statusOfProblemInPrimal(int &lastCleaned, int type,
827  ClpSimplexProgress *progress,
828  bool doFactorization,
829  int ifValuesPass,
830  ClpSimplex *originalModel)
831{
832  int dummy; // for use in generalExpanded
833  int saveFirstFree = firstFree_;
834  double saveOriginalWeight = infeasibilityCost_;
835  // number of pivots done
836  int numberPivots = factorization_->pivots();
837  if (type == 2) {
838    // trouble - restore solution
839    CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
840    CoinMemcpyN(savedSolution_ + numberColumns_,
841      numberRows_, rowActivityWork_);
842    CoinMemcpyN(savedSolution_,
843      numberColumns_, columnActivityWork_);
844    // restore extra stuff
845    matrix_->generalExpanded(this, 6, dummy);
846    forceFactorization_ = 1; // a bit drastic but ..
847    pivotRow_ = -1; // say no weights update
849  }
850  int saveThreshold = factorization_->sparseThreshold();
851  int tentativeStatus = problemStatus_;
852  int numberThrownOut = 1; // to loop round on bad factorization in values pass
853  double lastSumInfeasibility = COIN_DBL_MAX;
854  if (numberIterations_)
855    lastSumInfeasibility = nonLinearCost_->sumInfeasibilities();
856  int nPass = 0;
857  while (numberThrownOut) {
858    int nSlackBasic = 0;
859    if (nPass) {
860      for (int i = 0; i < numberRows_; i++) {
861        if (getRowStatus(i) == basic)
862          nSlackBasic++;
863      }
864    }
865    nPass++;
866    if (problemStatus_ > -3 || problemStatus_ == -4) {
867      // factorize
868      // later on we will need to recover from singularities
869      // also we could skip if first time
870      // do weights
871      // This may save pivotRow_ for use
872      if (doFactorization)
873        primalColumnPivot_->saveWeights(this, 1);
874
875      if ((type && doFactorization) || nSlackBasic == numberRows_) {
876        // is factorization okay?
877        int factorStatus = internalFactorize(1);
878        if (factorStatus) {
879          if (solveType_ == 2 + 8) {
880            // say odd
881            problemStatus_ = 5;
882            return;
883          }
884          if (type != 1 || largestPrimalError_ > 1.0e3
885            || largestDualError_ > 1.0e3) {
886            // switch off dense
887            int saveDense = factorization_->denseThreshold();
888            factorization_->setDenseThreshold(0);
889            // Go to safe
890            factorization_->pivotTolerance(0.99);
891            // make sure will do safe factorization
892            pivotVariable_[0] = -1;
893            internalFactorize(2);
894            factorization_->setDenseThreshold(saveDense);
895            // restore extra stuff
896            matrix_->generalExpanded(this, 6, dummy);
897          } else {
898            // no - restore previous basis
899            // Keep any flagged variables
900            int i;
901            for (i = 0; i < numberRows_ + numberColumns_; i++) {
902              if (flagged(i))
903                saveStatus_[i] |= 64; //say flagged
904            }
905            CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
906            if (numberPivots <= 1) {
907              // throw out something
908              if (sequenceIn_ >= 0 && getStatus(sequenceIn_) != basic) {
909                setFlagged(sequenceIn_);
910              } else if (sequenceOut_ >= 0 && getStatus(sequenceOut_) != basic) {
911                setFlagged(sequenceOut_);
912              }
913              double newTolerance = CoinMax(0.5 + 0.499 * randomNumberGenerator_.randomDouble(), factorization_->pivotTolerance());
914              factorization_->pivotTolerance(newTolerance);
915            } else {
916              // Go to safe
917              factorization_->pivotTolerance(0.99);
918            }
919            CoinMemcpyN(savedSolution_ + numberColumns_,
920              numberRows_, rowActivityWork_);
921            CoinMemcpyN(savedSolution_,
922              numberColumns_, columnActivityWork_);
923            // restore extra stuff
924            matrix_->generalExpanded(this, 6, dummy);
925            matrix_->generalExpanded(this, 5, dummy);
926            forceFactorization_ = 1; // a bit drastic but ..
927            type = 2;
928            if (internalFactorize(2) != 0) {
929              largestPrimalError_ = 1.0e4; // force other type
930            }
931          }
933        }
934      }
935      if (problemStatus_ != -4)
936        problemStatus_ = -3;
937    }
938    if (progress->realInfeasibility_[0] < 1.0e-1 && primalTolerance_ == 1.0e-7 && progress->iterationNumber_[0] > 0 && progress->iterationNumber_[CLP_PROGRESS - 1] - progress->iterationNumber_[0] > 25) {
939      // default - so user did not set
940      int iP;
941      double minAverage = COIN_DBL_MAX;
942      double maxAverage = 0.0;
943      for (iP = 0; iP < CLP_PROGRESS; iP++) {
944        int n = progress->numberInfeasibilities_[iP];
945        if (!n) {
946          break;
947        } else {
948          double average = progress->realInfeasibility_[iP];
949          if (average > 0.1)
950            break;
951          average /= static_cast< double >(n);
952          minAverage = CoinMin(minAverage, average);
953          maxAverage = CoinMax(maxAverage, average);
954        }
955      }
956      if (iP == CLP_PROGRESS && minAverage < 1.0e-5 && maxAverage < 1.0e-3) {
957        // change tolerance
958#if CBC_USEFUL_PRINTING > 0
959        printf("CCchanging tolerance\n");
960#endif
961        primalTolerance_ = 1.0e-6;
962        dblParam_[ClpPrimalTolerance] = 1.0e-6;
963        moreSpecialOptions_ |= 4194304;
964      }
965    }
966    // at this stage status is -3 or -5 if looks unbounded
967    // get primal and dual solutions
968    // put back original costs and then check
969    // createRim(4); // costs do not change
970    // May need to do more if column generation
971    dummy = 4;
972    matrix_->generalExpanded(this, 9, dummy);
973#ifndef CLP_CAUTION
974#define CLP_CAUTION 1
975#endif
976#if CLP_CAUTION
977    double lastAverageInfeasibility = sumDualInfeasibilities_ / static_cast< double >(numberDualInfeasibilities_ + 10);
978#endif
979#ifdef CLP_USER_DRIVEN
980    int status = eventHandler_->event(ClpEventHandler::goodFactorization);
981    if (status >= 0) {
982      lastSumInfeasibility = COIN_DBL_MAX;
983    }
984#endif
985    if (ifValuesPass && firstFree_ < 0) {
986      largestPrimalError_ = 1.0e7;
987      largestDualError_ = 1.0e7;
988    }
989    bool sayValuesPass = (firstFree_ >= 0);
990    //if (ifValuesPass&&numberIterations_==baseIteration_)
991    //sayValuesPass=true;
992    numberThrownOut = gutsOfSolution(NULL, NULL, sayValuesPass);
993    double sumInfeasibility = nonLinearCost_->sumInfeasibilities();
994    int reason2 = 0;
995#if CLP_CAUTION
996#if CLP_CAUTION == 2
997    double test2 = 1.0e5;
998#else
999    double test2 = 1.0e-1;
1000#endif
1001    if (!lastSumInfeasibility && sumInfeasibility && lastAverageInfeasibility < test2 && numberPivots > 10)
1002      reason2 = 3;
1003    if (lastSumInfeasibility < 1.0e-6 && sumInfeasibility > 1.0e-3 && numberPivots > 10)
1004      reason2 = 4;
1005#endif
1006    if (numberThrownOut)
1007      reason2 = 1;
1008    if ((sumInfeasibility > 1.0e7 && sumInfeasibility > 100.0 * lastSumInfeasibility
1009          && factorization_->pivotTolerance() < 0.11)
1010      || (largestPrimalError_ > 1.0e10 && largestDualError_ > 1.0e10))
1011      reason2 = 2;
1012    if (reason2) {
1013      problemStatus_ = tentativeStatus;
1014      doFactorization = true;
1015      if (numberPivots || (numberThrownOut == -123456789 && saveStatus_)) {
1016        // go back
1017        // trouble - restore solution
1018        CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
1019        CoinMemcpyN(savedSolution_ + numberColumns_,
1020          numberRows_, rowActivityWork_);
1021        CoinMemcpyN(savedSolution_,
1022          numberColumns_, columnActivityWork_);
1023        // restore extra stuff
1024        matrix_->generalExpanded(this, 6, dummy);
1025        if (reason2 < 3) {
1026          // Go to safe
1027          factorization_->pivotTolerance(CoinMin(0.99, 1.01 * factorization_->pivotTolerance()));
1028          forceFactorization_ = 1; // a bit drastic but ..
1029        } else if (forceFactorization_ < 0) {
1030          forceFactorization_ = CoinMin(numberPivots / 2, 100);
1031        } else {
1032          forceFactorization_ = CoinMin(forceFactorization_,
1033            CoinMax(3, numberPivots / 2));
1034        }
1035        pivotRow_ = -1; // say no weights update
1037        if (numberPivots == 1) {
1038          // throw out something
1039          if (sequenceIn_ >= 0 && getStatus(sequenceIn_) != basic) {
1040            setFlagged(sequenceIn_);
1041          } else if (sequenceOut_ >= 0 && getStatus(sequenceOut_) != basic) {
1042            setFlagged(sequenceOut_);
1043          }
1044        }
1045        type = 2; // so will restore weights
1046        if (internalFactorize(2) != 0) {
1047          largestPrimalError_ = 1.0e4; // force other type
1048        }
1049        numberPivots = 0;
1050        numberThrownOut = gutsOfSolution(NULL, NULL, (firstFree_ >= 0));
1051        //assert (!numberThrownOut);
1052#ifdef CLP_USEFUL_PRINTOUT
1053        if (numberThrownOut)
1054          printf("OUCH! - %d thrownout at %s line %d\n",
1055            numberThrownOut, __FILE__, __LINE__);
1056#endif
1057        sumInfeasibility = nonLinearCost_->sumInfeasibilities();
1058      }
1059    }
1060  }
1061  // Double check reduced costs if no action
1062  if (progress->lastIterationNumber(0) == numberIterations_) {
1063    if (primalColumnPivot_->looksOptimal()) {
1064      numberDualInfeasibilities_ = 0;
1065      sumDualInfeasibilities_ = 0.0;
1066    }
1067  }
1068  // If in primal and small dj give up
1069  if ((specialOptions_ & 1024) != 0 && !numberPrimalInfeasibilities_ && numberDualInfeasibilities_) {
1070    double average = sumDualInfeasibilities_ / (static_cast< double >(numberDualInfeasibilities_));
1071    if (numberIterations_ > 300 && average < 1.0e-4) {
1072      numberDualInfeasibilities_ = 0;
1073      sumDualInfeasibilities_ = 0.0;
1074    }
1075  }
1076  // Check if looping
1077  int loop;
1078  if (type != 2 && !ifValuesPass)
1079    loop = progress->looping();
1080  else
1081    loop = -1;
1082  if (loop >= 0) {
1083    if (!problemStatus_) {
1084      // declaring victory
1085      numberPrimalInfeasibilities_ = 0;
1086      sumPrimalInfeasibilities_ = 0.0;
1087    } else {
1088      problemStatus_ = loop; //exit if in loop
1089      problemStatus_ = 10; // instead - try other algorithm
1090      numberPrimalInfeasibilities_ = nonLinearCost_->numberInfeasibilities();
1091    }
1092    problemStatus_ = 10; // instead - try other algorithm
1093    return;
1094  } else if (loop < -1) {
1095    // Is it time for drastic measures
1096    if (nonLinearCost_->numberInfeasibilities() && progress->badTimes() > 5 && progress->oddState() < 10 && progress->oddState() >= 0) {
1097      progress->newOddState();
1098      nonLinearCost_->zapCosts();
1099    }
1100    // something may have changed
1101    gutsOfSolution(NULL, NULL, ifValuesPass != 0);
1102  }
1103  // If progress then reset costs
1104  if (loop == -1 && !nonLinearCost_->numberInfeasibilities() && progress->oddState() < 0) {
1105    createRim(4, false); // costs back
1106    delete nonLinearCost_;
1107    nonLinearCost_ = new ClpNonLinearCost(this);
1108    progress->endOddState();
1109    gutsOfSolution(NULL, NULL, ifValuesPass != 0);
1110  }
1111  // Flag to say whether to go to dual to clean up
1112  bool goToDual = false;
1114  //if((progressFlag_&2)!=0)
1115  //problemStatus_=-1;;
1116  progressFlag_ = 0; //reset progress flag
1117
1118  handler_->message(CLP_SIMPLEX_STATUS, messages_)
1119    << numberIterations_ << nonLinearCost_->feasibleReportCost();
1120  handler_->printing(nonLinearCost_->numberInfeasibilities() > 0)
1121    << nonLinearCost_->sumInfeasibilities() << nonLinearCost_->numberInfeasibilities();
1122  handler_->printing(sumDualInfeasibilities_ > 0.0)
1123    << sumDualInfeasibilities_ << numberDualInfeasibilities_;
1124  handler_->printing(numberDualInfeasibilitiesWithoutFree_
1125    < numberDualInfeasibilities_)
1126    << numberDualInfeasibilitiesWithoutFree_;
1127  handler_->message() << CoinMessageEol;
1128  if (!primalFeasible()) {
1129    nonLinearCost_->checkInfeasibilities(primalTolerance_);
1130    gutsOfSolution(NULL, NULL, ifValuesPass != 0);
1131    nonLinearCost_->checkInfeasibilities(primalTolerance_);
1132  }
1133  if (nonLinearCost_->numberInfeasibilities() > 0 && !progress->initialWeight_ && !ifValuesPass && infeasibilityCost_ == 1.0e10) {
1134    // first time infeasible - start up weight computation
1135    double *oldDj = dj_;
1136    double *oldCost = cost_;
1137    int numberRows2 = numberRows_ + numberExtraRows_;
1138    int numberTotal = numberRows2 + numberColumns_;
1139    dj_ = new double[numberTotal];
1140    cost_ = new double[numberTotal];
1141    reducedCostWork_ = dj_;
1142    rowReducedCost_ = dj_ + numberColumns_;
1143    objectiveWork_ = cost_;
1144    rowObjectiveWork_ = cost_ + numberColumns_;
1145    double direction = optimizationDirection_ * objectiveScale_;
1146    const double *obj = objective();
1147    memset(rowObjectiveWork_, 0, numberRows_ * sizeof(double));
1148    int iSequence;
1149    if (columnScale_)
1150      for (iSequence = 0; iSequence < numberColumns_; iSequence++)
1151        cost_[iSequence] = obj[iSequence] * direction * columnScale_[iSequence];
1152    else
1153      for (iSequence = 0; iSequence < numberColumns_; iSequence++)
1154        cost_[iSequence] = obj[iSequence] * direction;
1155    computeDuals(NULL);
1156    int numberSame = 0;
1157    int numberDifferent = 0;
1158    int numberZero = 0;
1159    int numberFreeSame = 0;
1160    int numberFreeDifferent = 0;
1161    int numberFreeZero = 0;
1162    int n = 0;
1163    for (iSequence = 0; iSequence < numberTotal; iSequence++) {
1164      if (getStatus(iSequence) != basic && !flagged(iSequence)) {
1165        // not basic
1166        double distanceUp = upper_[iSequence] - solution_[iSequence];
1167        double distanceDown = solution_[iSequence] - lower_[iSequence];
1168        double feasibleDj = dj_[iSequence];
1169        double infeasibleDj = oldDj[iSequence] - feasibleDj;
1170        double value = feasibleDj * infeasibleDj;
1171        if (distanceUp > primalTolerance_) {
1172          // Check if "free"
1173          if (distanceDown > primalTolerance_) {
1174            // free
1175            if (value > dualTolerance_) {
1176              numberFreeSame++;
1177            } else if (value < -dualTolerance_) {
1178              numberFreeDifferent++;
1179              dj_[n++] = feasibleDj / infeasibleDj;
1180            } else {
1181              numberFreeZero++;
1182            }
1183          } else {
1184            // should not be negative
1185            if (value > dualTolerance_) {
1186              numberSame++;
1187            } else if (value < -dualTolerance_) {
1188              numberDifferent++;
1189              dj_[n++] = feasibleDj / infeasibleDj;
1190            } else {
1191              numberZero++;
1192            }
1193          }
1194        } else if (distanceDown > primalTolerance_) {
1195          // should not be positive
1196          if (value > dualTolerance_) {
1197            numberSame++;
1198          } else if (value < -dualTolerance_) {
1199            numberDifferent++;
1200            dj_[n++] = feasibleDj / infeasibleDj;
1201          } else {
1202            numberZero++;
1203          }
1204        }
1205      }
1206      progress->initialWeight_ = -1.0;
1207    }
1208#ifdef CLP_USEFUL_PRINTOUT
1209    printf("XXXX %d same, %d different, %d zero, -- free %d %d %d\n",
1210      numberSame, numberDifferent, numberZero,
1211      numberFreeSame, numberFreeDifferent, numberFreeZero);
1212#endif
1213    // we want most to be same
1214    if (n) {
1215      std::sort(dj_, dj_ + n);
1216#ifdef CLP_USEFUL_PRINTOUT
1217      double most = 0.95;
1218      int which = static_cast< int >((1.0 - most) * static_cast< double >(n));
1219      double take2 = -dj_[which] * infeasibilityCost_;
1220      printf("XXXX inf cost %g take %g (range %g %g)\n", infeasibilityCost_, take2, -dj_[0] * infeasibilityCost_, -dj_[n - 1] * infeasibilityCost_);
1221#endif
1222      double take = -dj_[0] * infeasibilityCost_;
1223      // was infeasibilityCost_ = CoinMin(CoinMax(1000.0 * take, 1.0e8), 1.0000001e10);
1224      infeasibilityCost_ = CoinMin(CoinMax(1000.0 * take, 1.0e3), 1.0000001e10);
1225#ifdef CLP_USEFUL_PRINTOUT
1226      printf("XXXX changing weight to %g\n", infeasibilityCost_);
1227#endif
1228    }
1229    delete[] dj_;
1230    delete[] cost_;
1231    dj_ = oldDj;
1232    cost_ = oldCost;
1233    reducedCostWork_ = dj_;
1234    rowReducedCost_ = dj_ + numberColumns_;
1235    objectiveWork_ = cost_;
1236    rowObjectiveWork_ = cost_ + numberColumns_;
1237    if (n || matrix_->type() >= 15)
1238      gutsOfSolution(NULL, NULL, ifValuesPass != 0);
1239  }
1240  double trueInfeasibility = nonLinearCost_->sumInfeasibilities();
1241  if (!nonLinearCost_->numberInfeasibilities() && infeasibilityCost_ == 1.0e10 && !ifValuesPass && true) {
1242    // relax if default
1243    infeasibilityCost_ = CoinMin(CoinMax(100.0 * sumDualInfeasibilities_, 1.0e8), 1.00000001e10);
1244    // reset looping criterion
1245    progress->reset();
1246    trueInfeasibility = 1.123456e10;
1247  }
1248  if (trueInfeasibility > 1.0) {
1249    // If infeasibility going up may change weights
1250    double testValue = trueInfeasibility - 1.0e-4 * (10.0 + trueInfeasibility);
1251    double lastInf = progress->lastInfeasibility(1);
1252    double lastInf3 = progress->lastInfeasibility(3);
1253    double thisObj = progress->lastObjective(0);
1254    double thisInf = progress->lastInfeasibility(0);
1255    thisObj += infeasibilityCost_ * 2.0 * thisInf;
1256    double lastObj = progress->lastObjective(1);
1257    lastObj += infeasibilityCost_ * 2.0 * lastInf;
1258    double lastObj3 = progress->lastObjective(3);
1259    lastObj3 += infeasibilityCost_ * 2.0 * lastInf3;
1260    assert(thisObj <= COIN_DBL_MAX);
1261    assert(lastObj <= COIN_DBL_MAX);
1262    assert(lastObj3 <= COIN_DBL_MAX);
1263    if (lastObj < thisObj - 1.0e-5 * CoinMax(fabs(thisObj), fabs(lastObj)) - 1.0e-7
1264      && firstFree_ < 0 && thisInf >= lastInf) {
1265      if (handler_->logLevel() == 63)
1266        printf("lastobj %g this %g force %d\n", lastObj, thisObj, forceFactorization_);
1267      int maxFactor = factorization_->maximumPivots();
1268      if (maxFactor > 10) {
1269        if (forceFactorization_ < 0)
1270          forceFactorization_ = maxFactor;
1271        forceFactorization_ = CoinMax(1, (forceFactorization_ >> 2));
1272        if (handler_->logLevel() == 63)
1273          printf("Reducing factorization frequency to %d\n", forceFactorization_);
1274      }
1275    } else if (lastObj3 < thisObj - 1.0e-5 * CoinMax(fabs(thisObj), fabs(lastObj3)) - 1.0e-7
1276      && firstFree_ < 0 && thisInf >= lastInf) {
1277      if (handler_->logLevel() == 63)
1278        printf("lastobj3 %g this3 %g force %d\n", lastObj3, thisObj, forceFactorization_);
1279      int maxFactor = factorization_->maximumPivots();
1280      if (maxFactor > 10) {
1281        if (forceFactorization_ < 0)
1282          forceFactorization_ = maxFactor;
1283        forceFactorization_ = CoinMax(1, (forceFactorization_ * 2) / 3);
1284        if (handler_->logLevel() == 63)
1285          printf("Reducing factorization frequency to %d\n", forceFactorization_);
1286      }
1287    } else if (lastInf < testValue || trueInfeasibility == 1.123456e10) {
1288      if (infeasibilityCost_ < 1.0e14) {
1289        infeasibilityCost_ *= 1.5;
1290        // reset looping criterion
1291        progress->reset();
1292        if (handler_->logLevel() == 63)
1293          printf("increasing weight to %g\n", infeasibilityCost_);
1294        gutsOfSolution(NULL, NULL, ifValuesPass != 0);
1295      }
1296    }
1297  }
1298  // we may wish to say it is optimal even if infeasible
1299  bool alwaysOptimal = (specialOptions_ & 1) != 0;
1300  // give code benefit of doubt
1301  if (sumOfRelaxedDualInfeasibilities_ == 0.0 && sumOfRelaxedPrimalInfeasibilities_ == 0.0 && progress->objective_[CLP_PROGRESS - 1] >= progress->objective_[CLP_PROGRESS - 2] - 1.0e-9 * (10.0 + fabs(objectiveValue_))) {
1302    // say optimal (with these bounds etc)
1303    numberDualInfeasibilities_ = 0;
1304    sumDualInfeasibilities_ = 0.0;
1305    numberPrimalInfeasibilities_ = 0;
1306    sumPrimalInfeasibilities_ = 0.0;
1307    // But check if in sprint
1308    if (originalModel) {
1309      // Carry on and re-do
1310      numberDualInfeasibilities_ = -776;
1311    }
1312    // But if real primal infeasibilities nonzero carry on
1313    if (nonLinearCost_->numberInfeasibilities()) {
1314      // most likely to happen if infeasible
1315      double relaxedToleranceP = primalTolerance_;
1316      // we can't really trust infeasibilities if there is primal error
1317      double error = CoinMin(1.0e-2, largestPrimalError_);
1318      // allow tolerance at least slightly bigger than standard
1319      relaxedToleranceP = relaxedToleranceP + error;
1320      int ninfeas = nonLinearCost_->numberInfeasibilities();
1321      double sum = nonLinearCost_->sumInfeasibilities();
1322      double average = sum / static_cast< double >(ninfeas);
1323#ifdef COIN_DEVELOP
1324      if (handler_->logLevel() > 0)
1325        printf("nonLinearCost says infeasible %d summing to %g\n",
1326          ninfeas, sum);
1327#endif
1328      if (average > relaxedToleranceP) {
1329        sumOfRelaxedPrimalInfeasibilities_ = sum;
1330        numberPrimalInfeasibilities_ = ninfeas;
1331        sumPrimalInfeasibilities_ = sum;
1332#ifdef COIN_DEVELOP
1333        bool unflagged =
1334#endif
1335          unflag();
1336#ifdef COIN_DEVELOP
1337        if (unflagged && handler_->logLevel() > 0)
1338          printf(" - but flagged variables\n");
1339#endif
1340      }
1341    }
1342  }
1343  bool looksOptimal = (!numberDualInfeasibilities_ && !nonLinearCost_->sumInfeasibilities());
1344  // had ||(type==3&&problemStatus_!=-5) -- ??? why ????
1345  if ((dualFeasible() || problemStatus_ == -4) && (!ifValuesPass || looksOptimal || firstFree_ < 0)) {
1346    // see if extra helps
1347    if (nonLinearCost_->numberInfeasibilities() && (nonLinearCost_->sumInfeasibilities() > 1.0e-3 || sumOfRelaxedPrimalInfeasibilities_)
1348      /*&& alwaysOptimal*/) {
1349      //may need infeasiblity cost changed
1350      // we can see if we can construct a ray
1351      // make up a new objective
1352      double saveWeight = infeasibilityCost_;
1353      // save nonlinear cost as we are going to switch off costs
1354      //ClpNonLinearCost * nonLinear = nonLinearCost_;
1355      // do twice to make sure Primal solution has settled
1356      // put non-basics to bounds in case tolerance moved
1357      // put back original costs
1358      createRim(4);
1359      nonLinearCost_->checkInfeasibilities(0.0);
1360      gutsOfSolution(NULL, NULL, ifValuesPass != 0);
1361
1362      //infeasibilityCost_ = 1.0e100;
1363      // put back original costs
1364      createRim(4);
1365      nonLinearCost_->checkInfeasibilities(primalTolerance_);
1366      // may have fixed infeasibilities - double check
1367      if (nonLinearCost_->numberInfeasibilities() == 0) {
1368        // carry on
1369        problemStatus_ = -1;
1370        infeasibilityCost_ = saveWeight;
1371        nonLinearCost_->checkInfeasibilities(primalTolerance_);
1372      } else {
1373#ifndef MAX_INFEASIBILITY_COST
1374#define MAX_INFEASIBILITY_COST 1.0e18
1375#endif
1376        infeasibilityCost_ = 1.0e30;
1377        gutsOfSolution(NULL, NULL, ifValuesPass != 0 && firstFree_ >= 0);
1378        infeasibilityCost_ = CoinMax(saveWeight, saveOriginalWeight);
1379        if ((infeasibilityCost_ >= MAX_INFEASIBILITY_COST || numberDualInfeasibilities_ == 0) && perturbation_ == 101) {
1380          goToDual = unPerturb(); // stop any further perturbation
1381          if (nonLinearCost_->sumInfeasibilities() > 1.0e-1)
1382            goToDual = false;
1383          nonLinearCost_->checkInfeasibilities(primalTolerance_);
1384          numberDualInfeasibilities_ = 1; // carry on
1385          problemStatus_ = -1;
1386        } else if (numberDualInfeasibilities_ == 0 && largestDualError_ > 1.0e-2 && (moreSpecialOptions_ & (256 | 8192)) == 0) {
1387          goToDual = true;
1388          factorization_->pivotTolerance(CoinMax(0.9, factorization_->pivotTolerance()));
1389        }
1390        if (!goToDual) {
1391          if (infeasibilityCost_ >= MAX_INFEASIBILITY_COST || numberDualInfeasibilities_ == 0) {
1392            // we are infeasible - use as ray
1393            delete[] ray_;
1394            //if ((specialOptions_&(32|0x01000000))!=0x01000000) {
1395            if ((specialOptions_ & 0x03000000) == 0) {
1396              ray_ = new double[numberRows_];
1397              double *saveObjective = CoinCopyOfArray(objective(),
1398                numberColumns_);
1399              memset(objective(), 0, numberColumns_ * sizeof(double));
1400              infeasibilityCost_ = 1.0;
1401              createRim(4);
1402              nonLinearCost_->checkInfeasibilities(primalTolerance_);
1403              gutsOfSolution(NULL, NULL, false);
1404              memcpy(objective(), saveObjective, numberColumns_ * sizeof(double));
1405              delete[] saveObjective;
1406              // swap sign
1407              for (int i = 0; i < numberRows_; i++)
1408                ray_[i] = -dual_[i];
1409            } else {
1410              ray_ = NULL;
1411            }
1412#ifdef PRINT_RAY_METHOD
1413            printf("Primal creating infeasibility ray\n");
1414#endif
1415            // and get feasible duals
1416            infeasibilityCost_ = 0.0;
1417            createRim(4);
1418            nonLinearCost_->checkInfeasibilities(primalTolerance_);
1419            gutsOfSolution(NULL, NULL, ifValuesPass != 0);
1420            // so will exit
1421            infeasibilityCost_ = 1.0e30;
1422            // reset infeasibilities
1423            sumPrimalInfeasibilities_ = nonLinearCost_->sumInfeasibilities();
1424            ;
1425            numberPrimalInfeasibilities_ = nonLinearCost_->numberInfeasibilities();
1426          }
1427          if (infeasibilityCost_ < MAX_INFEASIBILITY_COST) {
1428            double testValue = MAX_INFEASIBILITY_COST;
1429            if (testValue == 1.0e18) {
1430              // Check it is not just noise
1431              const double *obj = objective();
1432              double largestCost = 0.0;
1433              if (columnScale_) {
1434                for (int i = 0; i < numberColumns_; i++) {
1435                  largestCost = CoinMax(largestCost, fabs(obj[i] * columnScale_[i]));
1436                }
1437              } else {
1438                for (int i = 0; i < numberColumns_; i++) {
1439                  largestCost = CoinMax(largestCost, fabs(obj[i]));
1440                }
1441              }
1442              testValue = 1.0e12 * (largestCost + 1.0e-6);
1443              if (numberDualInfeasibilities_) {
1444                double average = sumDualInfeasibilities_ / numberDualInfeasibilities_;
1445                testValue = CoinMax(testValue, average);
1446              }
1447              testValue = CoinMin(testValue, MAX_INFEASIBILITY_COST);
1448            }
1449            if (infeasibilityCost_ < testValue) {
1450              infeasibilityCost_ *= 5.0;
1451              // reset looping criterion
1452              progress->reset();
1454              handler_->message(CLP_PRIMAL_WEIGHT, messages_)
1455                << infeasibilityCost_
1456                << CoinMessageEol;
1457              // put back original costs and then check
1458              createRim(4);
1459              nonLinearCost_->checkInfeasibilities(0.0);
1460              gutsOfSolution(NULL, NULL, ifValuesPass != 0);
1461              problemStatus_ = -1; //continue
1462              goToDual = false;
1463            } else {
1464              // say infeasible
1465              problemStatus_ = 1;
1466            }
1467          } else {
1468            // say infeasible
1469            problemStatus_ = 1;
1470          }
1471        }
1472      }
1473    } else {
1474      // may be optimal
1475      if (perturbation_ == 101) {
1476        goToDual = unPerturb(); // stop any further perturbation
1477        if ((numberRows_ > 20000 || numberDualInfeasibilities_) && !numberTimesOptimal_)
1478          goToDual = false; // Better to carry on a bit longer
1479        lastCleaned = -1; // carry on
1480      }
1481      bool unflagged = (unflag() != 0);
1482      if (lastCleaned != numberIterations_ || unflagged) {
1483        handler_->message(CLP_PRIMAL_OPTIMAL, messages_)
1484          << primalTolerance_
1485          << CoinMessageEol;
1486        if (numberTimesOptimal_ < 4) {
1487          numberTimesOptimal_++;
1489          if (numberTimesOptimal_ == 1) {
1490            // better to have small tolerance even if slower
1491            factorization_->zeroTolerance(CoinMin(factorization_->zeroTolerance(), 1.0e-15));
1492          }
1493          lastCleaned = numberIterations_;
1494          if (primalTolerance_ != dblParam_[ClpPrimalTolerance])
1495            handler_->message(CLP_PRIMAL_ORIGINAL, messages_)
1496              << CoinMessageEol;
1497          double oldTolerance = primalTolerance_;
1498          primalTolerance_ = dblParam_[ClpPrimalTolerance];
1499#if 0
1500                         double * xcost = new double[numberRows_+numberColumns_];
1501                         double * xlower = new double[numberRows_+numberColumns_];
1502                         double * xupper = new double[numberRows_+numberColumns_];
1503                         double * xdj = new double[numberRows_+numberColumns_];
1504                         double * xsolution = new double[numberRows_+numberColumns_];
1505                         CoinMemcpyN(cost_, (numberRows_ + numberColumns_), xcost);
1506                         CoinMemcpyN(lower_, (numberRows_ + numberColumns_), xlower);
1507                         CoinMemcpyN(upper_, (numberRows_ + numberColumns_), xupper);
1508                         CoinMemcpyN(dj_, (numberRows_ + numberColumns_), xdj);
1509                         CoinMemcpyN(solution_, (numberRows_ + numberColumns_), xsolution);
1510#endif
1511          // put back original costs and then check
1512          createRim(4);
1513          nonLinearCost_->checkInfeasibilities(oldTolerance);
1514#if 0
1515                         int i;
1516                         for (i = 0; i < numberRows_ + numberColumns_; i++) {
1517                              if (cost_[i] != xcost[i])
1518                                   printf("** %d old cost %g new %g sol %g\n",
1519                                          i, xcost[i], cost_[i], solution_[i]);
1520                              if (lower_[i] != xlower[i])
1521                                   printf("** %d old lower %g new %g sol %g\n",
1522                                          i, xlower[i], lower_[i], solution_[i]);
1523                              if (upper_[i] != xupper[i])
1524                                   printf("** %d old upper %g new %g sol %g\n",
1525                                          i, xupper[i], upper_[i], solution_[i]);
1526                              if (dj_[i] != xdj[i])
1527                                   printf("** %d old dj %g new %g sol %g\n",
1528                                          i, xdj[i], dj_[i], solution_[i]);
1529                              if (solution_[i] != xsolution[i])
1530                                   printf("** %d old solution %g new %g sol %g\n",
1531                                          i, xsolution[i], solution_[i], solution_[i]);
1532                         }
1533                         delete [] xcost;
1534                         delete [] xupper;
1535                         delete [] xlower;
1536                         delete [] xdj;
1537                         delete [] xsolution;
1538#endif
1539          gutsOfSolution(NULL, NULL, ifValuesPass != 0);
1540          if (sumOfRelaxedDualInfeasibilities_ == 0.0 && sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
1541            // say optimal (with these bounds etc)
1542            numberDualInfeasibilities_ = 0;
1543            sumDualInfeasibilities_ = 0.0;
1544            numberPrimalInfeasibilities_ = 0;
1545            sumPrimalInfeasibilities_ = 0.0;
1546          }
1547          if (dualFeasible() && !nonLinearCost_->numberInfeasibilities() && lastCleaned >= 0)
1548            problemStatus_ = 0;
1549          else
1550            problemStatus_ = -1;
1551        } else {
1552          problemStatus_ = 0; // optimal
1553          if (lastCleaned < numberIterations_) {
1554            handler_->message(CLP_SIMPLEX_GIVINGUP, messages_)
1555              << CoinMessageEol;
1556          }
1557        }
1558      } else {
1559        /* Previous code here mostly works but
1560                    sumOfRelaxed is rubbish in primal
1561                 - so give benefit of doubt still */
1562        double error = CoinMin(1.0e-4, largestPrimalError_);
1563        // allow bigger tolerance than standard
1564        double saveTolerance = primalTolerance_;
1565        primalTolerance_ = 2.0 * primalTolerance_ + error;
1566        nonLinearCost_->checkInfeasibilities(primalTolerance_);
1567        double relaxedSum = nonLinearCost_->sumInfeasibilities();
1568        // back
1569        primalTolerance_ = saveTolerance;
1570        nonLinearCost_->checkInfeasibilities(primalTolerance_);
1571        if (alwaysOptimal || !relaxedSum)
1572          problemStatus_ = 0; // optimal
1573        else
1574          problemStatus_ = 1; // infeasible
1575      }
1576    }
1577  } else {
1578    // see if looks unbounded
1579    if (problemStatus_ == -5) {
1580      if (nonLinearCost_->numberInfeasibilities()) {
1581        if (infeasibilityCost_ >= MAX_INFEASIBILITY_COST && perturbation_ == 101) {
1582          // back off weight
1583          infeasibilityCost_ = 1.0e13;
1584          // reset looping criterion
1585          progress->reset();
1586          unPerturb(); // stop any further perturbation
1587        }
1588        //we need infeasiblity cost changed
1589        if (infeasibilityCost_ < MAX_INFEASIBILITY_COST) {
1590          infeasibilityCost_ *= 5.0;
1591          // reset looping criterion
1592          progress->reset();
1594          handler_->message(CLP_PRIMAL_WEIGHT, messages_)
1595            << infeasibilityCost_
1596            << CoinMessageEol;
1597          // put back original costs and then check
1598          createRim(4);
1599          gutsOfSolution(NULL, NULL, ifValuesPass != 0);
1600          problemStatus_ = -1; //continue
1601        } else {
1602          // say infeasible
1603          problemStatus_ = 1;
1604          // we are infeasible - use as ray
1605          delete[] ray_;
1606          ray_ = new double[numberRows_];
1607          CoinMemcpyN(dual_, numberRows_, ray_);
1608        }
1609      } else {
1610        // say unbounded
1611        problemStatus_ = 2;
1612      }
1613    } else {
1614      // carry on
1615      problemStatus_ = -1;
1616      if (type == 3 && !ifValuesPass) {
1617        //bool unflagged =
1618        unflag();
1619        if (sumDualInfeasibilities_ < 1.0e-3 || (sumDualInfeasibilities_ / static_cast< double >(numberDualInfeasibilities_)) < 1.0e-5 || progress->lastIterationNumber(0) == numberIterations_) {
1620          if (!numberPrimalInfeasibilities_) {
1621            if (numberTimesOptimal_ < 4) {
1622              numberTimesOptimal_++;
1624            } else {
1625              problemStatus_ = 0;
1626              secondaryStatus_ = 5;
1627            }
1628          }
1629        }
1630      }
1631    }
1632  }
1633  if (problemStatus_ == 0) {
1634    double objVal = (nonLinearCost_->feasibleCost()
1635      + objective_->nonlinearOffset());
1636    objVal /= (objectiveScale_ * rhsScale_);
1637    double tol = 1.0e-10 * CoinMax(fabs(objVal), fabs(objectiveValue_)) + 1.0e-8;
1638    if (fabs(objVal - objectiveValue_) > tol) {
1639#ifdef COIN_DEVELOP
1640      if (handler_->logLevel() > 0)
1641        printf("nonLinearCost has feasible obj of %g, objectiveValue_ is %g\n",
1642          objVal, objectiveValue_);
1643#endif
1644      objectiveValue_ = objVal;
1645    }
1646  }
1647  // save extra stuff
1648  matrix_->generalExpanded(this, 5, dummy);
1649  if (type == 0 || type == 1) {
1650    if (type != 1 || !saveStatus_) {
1651      // create save arrays
1652      delete[] saveStatus_;
1653      delete[] savedSolution_;
1654      saveStatus_ = new unsigned char[numberRows_ + numberColumns_];
1655      savedSolution_ = new double[numberRows_ + numberColumns_];
1656    }
1657    // save arrays
1658    CoinMemcpyN(status_, numberColumns_ + numberRows_, saveStatus_);
1659    CoinMemcpyN(rowActivityWork_,
1660      numberRows_, savedSolution_ + numberColumns_);
1661    CoinMemcpyN(columnActivityWork_, numberColumns_, savedSolution_);
1662  }
1663  // see if in Cbc etc
1664  bool inCbcOrOther = (specialOptions_ & 0x03000000) != 0;
1665  bool disaster = false;
1666  if (disasterArea_ && inCbcOrOther && disasterArea_->check()) {
1667    disasterArea_->saveInfo();
1668    disaster = true;
1669  }
1670  if (disaster)
1671    problemStatus_ = 3;
1672  if (problemStatus_ < 0 && !changeMade_) {
1673    problemStatus_ = 4; // unknown
1674  }
1675  lastGoodIteration_ = numberIterations_;
1676  if (numberIterations_ > lastBadIteration_ + 100)
1677    moreSpecialOptions_ &= ~16; // clear check accuracy flag
1678  if ((moreSpecialOptions_ & 256) != 0)
1679    goToDual = false;
1680  if (goToDual || (numberIterations_ > 1000 && largestPrimalError_ > 1.0e6 && largestDualError_ > 1.0e6)) {
1681    problemStatus_ = 10; // try dual
1682    // See if second call
1683    if ((moreSpecialOptions_ & 256) != 0 && nonLinearCost_->sumInfeasibilities() > 1.0e2) {
1684      numberPrimalInfeasibilities_ = nonLinearCost_->numberInfeasibilities();
1685      sumPrimalInfeasibilities_ = nonLinearCost_->sumInfeasibilities();
1686      // say infeasible
1687      if (numberPrimalInfeasibilities_ && largestPrimalError_ < 1.0e-1)
1688        problemStatus_ = 1;
1689    }
1690  }
1691  // make sure first free monotonic
1692  if (firstFree_ >= 0 && saveFirstFree >= 0) {
1693    if (numberIterations_) {
1694      firstFree_ = saveFirstFree;
1695    } else {
1696      firstFree_ = -1;
1697      nextSuperBasic(1, NULL);
1698    }
1699  }
1700  if (doFactorization) {
1701    // restore weights (if saved) - also recompute infeasibility list
1702    if (tentativeStatus > -3)
1703      primalColumnPivot_->saveWeights(this, (type < 2) ? 2 : 4);
1704    else
1705      primalColumnPivot_->saveWeights(this, 3);
1706    if (saveThreshold) {
1707      // use default at present
1708      factorization_->sparseThreshold(0);
1709      factorization_->goSparse();
1710    }
1711  }
1712  if (problemStatus_ == 1) {
1713    // compute true objective value
1714    computeObjectiveValue(true);
1715  }
1716  // Allow matrices to be sorted etc
1717  int fake = -999; // signal sort
1718  matrix_->correctSequence(this, fake, fake);
1719}
1720/*
1721   Row array has pivot column
1722   This chooses pivot row.
1723   For speed, we may need to go to a bucket approach when many
1724   variables go through bounds
1725   On exit rhsArray will have changes in costs of basic variables
1726*/
1727void ClpSimplexPrimal::primalRow(CoinIndexedVector *rowArray,
1728  CoinIndexedVector *rhsArray,
1729  CoinIndexedVector *spareArray,
1730  int valuesPass)
1731{
1732  double saveDj = dualIn_;
1733  if (valuesPass && objective_->type() < 2) {
1734    dualIn_ = cost_[sequenceIn_];
1735
1736    double *work = rowArray->denseVector();
1737    int number = rowArray->getNumElements();
1738    int *which = rowArray->getIndices();
1739
1740    int iIndex;
1741    for (iIndex = 0; iIndex < number; iIndex++) {
1742
1743      int iRow = which[iIndex];
1744      double alpha = work[iIndex];
1745      int iPivot = pivotVariable_[iRow];
1746      dualIn_ -= alpha * cost_[iPivot];
1747    }
1748    // determine direction here
1749    if (dualIn_ < -dualTolerance_) {
1750      directionIn_ = 1;
1751    } else if (dualIn_ > dualTolerance_) {
1752      directionIn_ = -1;
1753    } else {
1754      // towards nearest bound
1755      if (valueIn_ - lowerIn_ < upperIn_ - valueIn_) {
1756        directionIn_ = -1;
1757        dualIn_ = dualTolerance_;
1758      } else {
1759        directionIn_ = 1;
1760        dualIn_ = -dualTolerance_;
1761      }
1762    }
1763  }
1764
1765  // sequence stays as row number until end
1766  pivotRow_ = -1;
1767  int numberRemaining = 0;
1768
1769  double totalThru = 0.0; // for when variables flip
1770  // Allow first few iterations to take tiny
1771  double acceptablePivot = 1.0e-1 * acceptablePivot_;
1772  if (numberIterations_ > 100)
1773    acceptablePivot = acceptablePivot_;
1774  if (factorization_->pivots() > 10)
1775    acceptablePivot = 1.0e+3 * acceptablePivot_; // if we have iterated be more strict
1776  else if (factorization_->pivots() > 5)
1777    acceptablePivot = 1.0e+2 * acceptablePivot_; // if we have iterated be slightly more strict
1778  else if (factorization_->pivots())
1779    acceptablePivot = acceptablePivot_; // relax
1780  double bestEverPivot = acceptablePivot;
1781  int lastPivotRow = -1;
1782  double lastPivot = 0.0;
1783  double lastTheta = 1.0e50;
1784
1785  // use spareArrays to put ones looked at in
1786  // First one is list of candidates
1787  // We could compress if we really know we won't need any more
1788  // Second array has current set of pivot candidates
1789  // with a backup list saved in double * part of indexed vector
1790
1791  // pivot elements
1792  double *spare;
1793  // indices
1794  int *index;
1795  spareArray->clear();
1796  spare = spareArray->denseVector();
1797  index = spareArray->getIndices();
1798
1799  // we also need somewhere for effective rhs
1800  double *rhs = rhsArray->denseVector();
1801  // and we can use indices to point to alpha
1802  // that way we can store fabs(alpha)
1803  int *indexPoint = rhsArray->getIndices();
1804  //int numberFlip=0; // Those which may change if flips
1805
1806  /*
1807       First we get a list of possible pivots.  We can also see if the
1808       problem looks unbounded.
1809
1810       At first we increase theta and see what happens.  We start
1811       theta at a reasonable guess.  If in right area then we do bit by bit.
1812       We save possible pivot candidates
1813
1814      */
1815
1816  // do first pass to get possibles
1817  // We can also see if unbounded
1818
1819  double *work = rowArray->denseVector();
1820  int number = rowArray->getNumElements();
1821  int *which = rowArray->getIndices();
1822
1823  // we need to swap sign if coming in from ub
1824  double way = directionIn_;
1825  double maximumMovement;
1826  if (way > 0.0)
1827    maximumMovement = CoinMin(1.0e30, upperIn_ - valueIn_);
1828  else
1829    maximumMovement = CoinMin(1.0e30, valueIn_ - lowerIn_);
1830
1831  double averageTheta = nonLinearCost_->averageTheta();
1832  double tentativeTheta = CoinMin(10.0 * averageTheta, maximumMovement);
1833  double upperTheta = maximumMovement;
1834  if (tentativeTheta > 0.5 * maximumMovement)
1835    tentativeTheta = maximumMovement;
1836  bool thetaAtMaximum = tentativeTheta == maximumMovement;
1837  // In case tiny bounds increase
1838  if (maximumMovement < 1.0)
1839    tentativeTheta *= 1.1;
1840  double dualCheck = fabs(dualIn_);
1841  // but make a bit more pessimistic
1842  dualCheck = CoinMax(dualCheck - 100.0 * dualTolerance_, 0.99 * dualCheck);
1843
1844  int iIndex;
1845  int pivotOne = -1;
1846  //#define CLP_DEBUG
1847#ifdef CLP_DEBUG
1848  if (numberIterations_ == -3839 || numberIterations_ == -3840) {
1849    double dj = cost_[sequenceIn_];
1850    printf("cost in on %d is %g, dual in %g\n", sequenceIn_, dj, dualIn_);
1851    for (iIndex = 0; iIndex < number; iIndex++) {
1852
1853      int iRow = which[iIndex];
1854      double alpha = work[iIndex];
1855      int iPivot = pivotVariable_[iRow];
1856      dj -= alpha * cost_[iPivot];
1857      printf("row %d var %d current %g %g %g, alpha %g so sol => %g (cost %g, dj %g)\n",
1858        iRow, iPivot, lower_[iPivot], solution_[iPivot], upper_[iPivot],
1859        alpha, solution_[iPivot] - 1.0e9 * alpha, cost_[iPivot], dj);
1860    }
1861  }
1862#endif
1863  while (true) {
1864    pivotOne = -1;
1865    totalThru = 0.0;
1866    //double totalThruFake=0.0;
1867    //int nThru=0;
1868    // We also re-compute reduced cost
1869    numberRemaining = 0;
1870    dualIn_ = cost_[sequenceIn_];
1871#ifndef NDEBUG
1872    //double tolerance = primalTolerance_ * 1.002;
1873#endif
1874    for (iIndex = 0; iIndex < number; iIndex++) {
1875
1876      int iRow = which[iIndex];
1877      double alpha = work[iIndex];
1878      int iPivot = pivotVariable_[iRow];
1879      // faster without if (cost_[iPivot])
1880      dualIn_ -= alpha * cost_[iPivot];
1881      alpha *= way;
1882      double oldValue = solution_[iPivot];
1883      // get where in bound sequence
1884      // note that after this alpha is actually fabs(alpha)
1885      bool possible;
1886      // do computation same way as later on in primal
1887      if (alpha > 0.0) {
1888        // basic variable going towards lower bound
1889        double bound = lower_[iPivot];
1890        // must be exactly same as when used
1891        double change = tentativeTheta * alpha;
1892        possible = (oldValue - change) <= bound + primalTolerance_;
1893        oldValue -= bound;
1894      } else {
1895        // basic variable going towards upper bound
1896        double bound = upper_[iPivot];
1897        // must be exactly same as when used
1898        double change = tentativeTheta * alpha;
1899        possible = (oldValue - change) >= bound - primalTolerance_;
1900        oldValue = bound - oldValue;
1901        alpha = -alpha;
1902      }
1903      double value;
1904      //assert (oldValue >= -10.0*tolerance);
1905      if (possible) {
1906        value = oldValue - upperTheta * alpha;
1907#ifdef CLP_USER_DRIVEN1
1908        if (!userChoiceValid1(this, iPivot, oldValue,
1909              upperTheta, alpha, work[iIndex] * way))
1910          value = 0.0; // say can't use
1911#endif
1912        if (value < -primalTolerance_ && alpha >= acceptablePivot) {
1913          upperTheta = (oldValue + primalTolerance_) / alpha;
1914          pivotOne = numberRemaining;
1915        }
1917        spare[numberRemaining] = alpha;
1918        rhs[numberRemaining] = oldValue;
1919        indexPoint[numberRemaining] = iIndex;
1920        index[numberRemaining++] = iRow;
1921        totalThru += alpha; // need more if dubious pricing
1922#if 0 // wasdef CLP_USER_DRIVEN
1923                    clpUserStruct info;
1924                    info.type=1;
1925                    info.row=iRow;
1926                    info.sequence=pivotVariable_[iRow];
1927                    info.alpha=way*work[iIndex];
1928                    double tempThru=totalThruFake+alpha;
1929                    info.totalThru=totalThruFake+alpha;
1930                    info.rhs=oldValue;
1931                    info.printing=printing;
1932                    eventHandler_->eventWithInfo(ClpEventHandler::pivotRow,
1933                                                 &info);
1934                    totalThruFake=info.totalThru;
1935                    if (fabs(totalThruFake-tempThru)>1.0e-4)
1936                      printf("%d thru - fake %g real %g\n",nThru,totalThruFake,totalThru);
1937                    nThru++;
1938#endif
1939        setActive(iRow);
1940        //} else if (value<primalTolerance_*1.002) {
1941        // May change if is a flip
1942        //indexRhs[numberFlip++]=iRow;
1943      }
1944    }
1945#if 0 // was def CLP_USER_DRIVEN
1946          clpUserStruct info;
1947          info.type=4;
1948          info.alpha=(upperTheta+1.0e-5)*way;
1949          eventHandler_->eventWithInfo(ClpEventHandler::pivotRow,
1950                                       &info);
1951          if (upperTheta < maximumMovement && totalThru*infeasibilityCost_ >= 1.0001 * dualCheck) {
1952            if (totalThruFake*infeasibilityCost_ < 1.0001 * dualCheck) {
1953              printf("fake allows through - it %d - fake %g real %g - dj %g\n",
1954                     numberIterations_,totalThruFake,totalThru,dualCheck);
1955            } else if (totalThruFake!=totalThru) {
1956              printf("it %d - fake %g real %g - dj %g - ratio %g\n",
1957                     numberIterations_,totalThruFake,totalThru,dualCheck,
1958                     dualCheck/totalThruFake);
1959            }
1960          }
1961#endif
1962    if (upperTheta < maximumMovement && totalThru * infeasibilityCost_ >= 1.0001 * dualCheck) {
1963      // Can pivot here
1964      break;
1965    } else if (!thetaAtMaximum) {
1966      //printf("Going round with average theta of %g\n",averageTheta);
1967      tentativeTheta = maximumMovement;
1968      thetaAtMaximum = true; // seems to be odd compiler error
1969    } else {
1970      break;
1971    }
1972  }
1973  totalThru = 0.0;
1974
1975  theta_ = maximumMovement;
1976
1977  bool goBackOne = false;
1978  if (objective_->type() > 1)
1979    dualIn_ = saveDj;
1980
1981  //printf("%d remain out of %d\n",numberRemaining,number);
1982  int iTry = 0;
1983#define MAXTRY 1000
1984  if (numberRemaining && upperTheta < maximumMovement) {
1985    // First check if previously chosen one will work
1986    if (pivotOne >= 0 && 0) {
1987      double thruCost = infeasibilityCost_ * spare[pivotOne];
1988      if (thruCost >= 0.99 * fabs(dualIn_))
1989        COIN_DETAIL_PRINT(printf("Could pivot on %d as change %g dj %g\n",
1990          index[pivotOne], thruCost, dualIn_));
1991      double alpha = spare[pivotOne];
1992      double oldValue = rhs[pivotOne];
1993      theta_ = oldValue / alpha;
1994      pivotRow_ = pivotOne;
1995      // Stop loop
1996      iTry = MAXTRY;
1997    }
1998
1999    // first get ratio with tolerance
2000    for (; iTry < MAXTRY; iTry++) {
2001
2002      upperTheta = maximumMovement;
2003      int iBest = -1;
2004      for (iIndex = 0; iIndex < numberRemaining; iIndex++) {
2005
2006        double alpha = spare[iIndex];
2007        double oldValue = rhs[iIndex];
2008        double value = oldValue - upperTheta * alpha;
2009
2010#ifdef CLP_USER_DRIVEN1
2011        int sequenceOut = pivotVariable_[index[iIndex]];
2012        if (!userChoiceValid1(this, sequenceOut, oldValue,
2013              upperTheta, alpha, 0.0))
2014          value = 0.0; // say can't use
2015#endif
2016        if (value < -primalTolerance_ && alpha >= acceptablePivot) {
2017          upperTheta = (oldValue + primalTolerance_) / alpha;
2018          iBest = iIndex; // just in case weird numbers
2019        }
2020      }
2021      // now look at best in this lot
2022      // But also see how infeasible small pivots will make
2023      double sumInfeasibilities = 0.0;
2024      double bestPivot = acceptablePivot;
2025      pivotRow_ = -1;
2026      for (iIndex = 0; iIndex < numberRemaining; iIndex++) {
2027
2028        int iRow = index[iIndex];
2029        double alpha = spare[iIndex];
2030        double oldValue = rhs[iIndex];
2031        double value = oldValue - upperTheta * alpha;
2032
2033        if (value <= 0 || iBest == iIndex) {
2034          // how much would it cost to go thru and modify bound
2035          double trueAlpha = way * work[indexPoint[iIndex]];
2036          int jSequence = pivotVariable_[iRow];
2037          //if (printing)
2038          //jSequence += 1000000;
2039          totalThru += nonLinearCost_->changeInCost(jSequence, trueAlpha, rhs[iIndex]);
2040#if 0 //was def CLP_USER_DRIVEN
2041                         clpUserStruct info;
2042                         info.type=2;
2043                         info.row=iRow;
2044                         info.sequence=pivotVariable_[iRow];
2045                         info.alpha=trueAlpha;
2046                         info.totalThru=totalThru;
2047                         info.rhs=rhs[iIndex];
2048                         info.printing=printing;
2049                         eventHandler_->eventWithInfo(ClpEventHandler::pivotRow,
2050                                                      &info);
2051                         totalThru=info.totalThru;
2052                         rhs[iIndex]=info.rhs;
2053#endif
2054          setActive(iRow);
2055          if (alpha > bestPivot) {
2056            bestPivot = alpha;
2057            theta_ = oldValue / bestPivot;
2058            pivotRow_ = iIndex;
2059          } else if (alpha < acceptablePivot
2060#ifdef CLP_USER_DRIVEN1
2061            || !userChoiceValid1(this, pivotVariable_[index[iIndex]],
2062                 oldValue, upperTheta, alpha, 0.0)
2063#endif
2064          ) {
2065            if (value < -primalTolerance_)
2066              sumInfeasibilities += -value - primalTolerance_;
2067          }
2068        }
2069      }
2070#if 0 // was def CLP_USER_DRIVEN
2071               clpUserStruct info;
2072               info.type=5;
2073               info.alpha=(theta_+1.0e-5)*way;
2074               eventHandler_->eventWithInfo(ClpEventHandler::pivotRow,
2075                                       &info);
2076#endif
2077      if (bestPivot < 0.1 * bestEverPivot && bestEverPivot > 1.0e-6 && bestPivot < 1.0e-3) {
2078        // back to previous one
2079        goBackOne = true;
2080        break;
2081      } else if (pivotRow_ == -1 && upperTheta > largeValue_) {
2082        if (lastPivot > acceptablePivot) {
2083          // back to previous one
2084          goBackOne = true;
2085        } else {
2086          // can only get here if all pivots so far too small
2087        }
2088        break;
2089      } else if (totalThru >= dualCheck) {
2090        if (sumInfeasibilities > primalTolerance_ && !nonLinearCost_->numberInfeasibilities()) {
2091          // Looks a bad choice
2092          if (lastPivot > acceptablePivot) {
2093            goBackOne = true;
2094          } else {
2095            // say no good
2096            dualIn_ = 0.0;
2097          }
2098        }
2099        break; // no point trying another loop
2100      } else {
2101        lastPivotRow = pivotRow_;
2102        lastTheta = theta_;
2103        if (bestPivot > bestEverPivot)
2104          bestEverPivot = bestPivot;
2105      }
2106    }
2107    // can get here without pivotRow_ set but with lastPivotRow
2108    if (goBackOne || (pivotRow_ < 0 && lastPivotRow >= 0)) {
2109      // back to previous one
2110      pivotRow_ = lastPivotRow;
2111      theta_ = lastTheta;
2112    }
2113  } else if (pivotRow_ < 0 && maximumMovement > 1.0e20) {
2114    // looks unbounded
2115    valueOut_ = COIN_DBL_MAX; // say odd
2116    if (nonLinearCost_->numberInfeasibilities()) {
2117      // but infeasible??
2118      // move variable but don't pivot
2119      tentativeTheta = 1.0e50;
2120      for (iIndex = 0; iIndex < number; iIndex++) {
2121        int iRow = which[iIndex];
2122        double alpha = work[iIndex];
2123        int iPivot = pivotVariable_[iRow];
2124        alpha *= way;
2125        double oldValue = solution_[iPivot];
2126        // get where in bound sequence
2127        // note that after this alpha is actually fabs(alpha)
2128        if (alpha > 0.0) {
2129          // basic variable going towards lower bound
2130          double bound = lower_[iPivot];
2131          oldValue -= bound;
2132        } else {
2133          // basic variable going towards upper bound
2134          double bound = upper_[iPivot];
2135          oldValue = bound - oldValue;
2136          alpha = -alpha;
2137        }
2138        if (oldValue - tentativeTheta * alpha < 0.0) {
2139          tentativeTheta = oldValue / alpha;
2140        }
2141      }
2142      // If free in then see if we can get to 0.0
2143      if (lowerIn_ < -1.0e20 && upperIn_ > 1.0e20) {
2144        if (dualIn_ * valueIn_ > 0.0) {
2145          if (fabs(valueIn_) < 1.0e-2 && (tentativeTheta < fabs(valueIn_) || tentativeTheta > 1.0e20)) {
2146            tentativeTheta = fabs(valueIn_);
2147          }
2148        }
2149      }
2150      if (tentativeTheta < 1.0e10)
2151        valueOut_ = valueIn_ + way * tentativeTheta;
2152    }
2153  }
2154  //if (iTry>50)
2155  //printf("** %d tries\n",iTry);
2156  if (pivotRow_ >= 0) {
2157    int position = pivotRow_; // position in list
2158    pivotRow_ = index[position];
2159    alpha_ = work[indexPoint[position]];
2160    // translate to sequence
2161    sequenceOut_ = pivotVariable_[pivotRow_];
2162    valueOut_ = solution(sequenceOut_);
2163    lowerOut_ = lower_[sequenceOut_];
2164    upperOut_ = upper_[sequenceOut_];
2165#define MINIMUMTHETA 1.0e-12
2166    // Movement should be minimum for anti-degeneracy - unless
2167    // fixed variable out
2168    double minimumTheta;
2169    if (upperOut_ > lowerOut_)
2170      minimumTheta = MINIMUMTHETA;
2171    else
2172      minimumTheta = 0.0;
2173    // But can't go infeasible
2174    double distance;
2175    if (alpha_ * way > 0.0)
2176      distance = valueOut_ - lowerOut_;
2177    else
2178      distance = upperOut_ - valueOut_;
2179    if (distance - minimumTheta * fabs(alpha_) < -primalTolerance_)
2180      minimumTheta = CoinMax(0.0, (distance + 0.5 * primalTolerance_) / fabs(alpha_));
2181    // will we need to increase tolerance
2182    //#define CLP_DEBUG
2183    double largestInfeasibility = primalTolerance_;
2184    if (theta_ < minimumTheta && (specialOptions_ & 4) == 0 && !valuesPass) {
2185      theta_ = minimumTheta;
2186      for (iIndex = 0; iIndex < numberRemaining - numberRemaining; iIndex++) {
2187        largestInfeasibility = CoinMax(largestInfeasibility,
2188          -(rhs[iIndex] - spare[iIndex] * theta_));
2189      }
2190//#define CLP_DEBUG
2191#ifdef CLP_DEBUG
2192      if (largestInfeasibility > primalTolerance_ && (handler_->logLevel() & 32) > -1)
2193        printf("Primal tolerance increased from %g to %g\n",
2194          primalTolerance_, largestInfeasibility);
2195#endif
2196      //#undef CLP_DEBUG
2197      primalTolerance_ = CoinMax(primalTolerance_, largestInfeasibility);
2198    }
2199    // Need to look at all in some cases
2200    if (theta_ > tentativeTheta) {
2201      for (iIndex = 0; iIndex < number; iIndex++)
2202        setActive(which[iIndex]);
2203    }
2204    if (way < 0.0)
2205      theta_ = -theta_;
2206    double newValue = valueOut_ - theta_ * alpha_;
2207    // If  4 bit set - Force outgoing variables to exact bound (primal)
2208    if (alpha_ * way < 0.0) {
2209      directionOut_ = -1; // to upper bound
2210      if (fabs(theta_) > 1.0e-6 || (specialOptions_ & 4) != 0) {
2211        upperOut_ = nonLinearCost_->nearest(sequenceOut_, newValue);
2212      } else {
2213        upperOut_ = newValue;
2214      }
2215    } else {
2216      directionOut_ = 1; // to lower bound
2217      if (fabs(theta_) > 1.0e-6 || (specialOptions_ & 4) != 0) {
2218        lowerOut_ = nonLinearCost_->nearest(sequenceOut_, newValue);
2219      } else {
2220        lowerOut_ = newValue;
2221      }
2222    }
2223    dualOut_ = reducedCost(sequenceOut_);
2224  } else if (maximumMovement < 1.0e20) {
2225    // flip
2226    pivotRow_ = -2; // so we can tell its a flip
2227    sequenceOut_ = sequenceIn_;
2228    valueOut_ = valueIn_;
2229    dualOut_ = dualIn_;
2230    lowerOut_ = lowerIn_;
2231    upperOut_ = upperIn_;
2232    alpha_ = 0.0;
2233    if (way < 0.0) {
2234      directionOut_ = 1; // to lower bound
2235      theta_ = lowerOut_ - valueOut_;
2236    } else {
2237      directionOut_ = -1; // to upper bound
2238      theta_ = upperOut_ - valueOut_;
2239    }
2240  }
2241
2242  double theta1 = CoinMax(theta_, 1.0e-12);
2243  double theta2 = numberIterations_ * nonLinearCost_->averageTheta();
2244  // Set average theta
2245  nonLinearCost_->setAverageTheta((theta1 + theta2) / (static_cast< double >(numberIterations_ + 1)));
2246  //if (numberIterations_%1000==0)
2247  //printf("average theta is %g\n",nonLinearCost_->averageTheta());
2248
2249  // clear arrays
2250
2251  CoinZeroN(spare, numberRemaining);
2252
2253  // put back original bounds etc
2254  CoinMemcpyN(index, numberRemaining, rhsArray->getIndices());
2255  rhsArray->setNumElements(numberRemaining);
2256  rhsArray->setPacked();
2257  nonLinearCost_->goBackAll(rhsArray);
2258  rhsArray->clear();
2259}
2260/*
2261   Chooses primal pivot column
2262   updateArray has cost updates (also use pivotRow_ from last iteration)
2263   Would be faster with separate region to scan
2264   and will have this (with square of infeasibility) when steepest
2265   For easy problems we can just choose one of the first columns we look at
2266*/
2268  CoinIndexedVector *spareRow1,
2269  CoinIndexedVector *spareRow2,
2270  CoinIndexedVector *spareColumn1,
2271  CoinIndexedVector *spareColumn2)
2272{
2273
2274  ClpMatrixBase *saveMatrix = matrix_;
2275  double *saveRowScale = rowScale_;
2276  if (scaledMatrix_) {
2277    rowScale_ = NULL;
2278    matrix_ = scaledMatrix_;
2279  }
2281    spareRow2, spareColumn1,
2282    spareColumn2);
2283  if (scaledMatrix_) {
2284    matrix_ = saveMatrix;
2285    rowScale_ = saveRowScale;
2286  }
2287  if (sequenceIn_ >= 0) {
2288    valueIn_ = solution_[sequenceIn_];
2289    dualIn_ = dj_[sequenceIn_];
2290    if (nonLinearCost_->lookBothWays()) {
2291      // double check
2292      ClpSimplex::Status status = getStatus(sequenceIn_);
2293
2294      switch (status) {
2295      case ClpSimplex::atUpperBound:
2296        if (dualIn_ < 0.0) {
2297          // move to other side
2298          COIN_DETAIL_PRINT(printf("For %d U (%g, %g, %g) dj changed from %g",
2299            sequenceIn_, lower_[sequenceIn_], solution_[sequenceIn_],
2300            upper_[sequenceIn_], dualIn_));
2301          dualIn_ -= nonLinearCost_->changeUpInCost(sequenceIn_);
2302          COIN_DETAIL_PRINT(printf(" to %g\n", dualIn_));
2303          nonLinearCost_->setOne(sequenceIn_, upper_[sequenceIn_] + 2.0 * currentPrimalTolerance());
2304          setStatus(sequenceIn_, ClpSimplex::atLowerBound);
2305        }
2306        break;
2307      case ClpSimplex::atLowerBound:
2308        if (dualIn_ > 0.0) {
2309          // move to other side
2310          COIN_DETAIL_PRINT(printf("For %d L (%g, %g, %g) dj changed from %g",
2311            sequenceIn_, lower_[sequenceIn_], solution_[sequenceIn_],
2312            upper_[sequenceIn_], dualIn_));
2313          dualIn_ -= nonLinearCost_->changeDownInCost(sequenceIn_);
2314          COIN_DETAIL_PRINT(printf(" to %g\n", dualIn_));
2315          nonLinearCost_->setOne(sequenceIn_, lower_[sequenceIn_] - 2.0 * currentPrimalTolerance());
2316          setStatus(sequenceIn_, ClpSimplex::atUpperBound);
2317        }
2318        break;
2319      default:
2320        break;
2321      }
2322    }
2323    lowerIn_ = lower_[sequenceIn_];
2324    upperIn_ = upper_[sequenceIn_];
2325    if (dualIn_ > 0.0)
2326      directionIn_ = -1;
2327    else
2328      directionIn_ = 1;
2329  } else {
2330    sequenceIn_ = -1;
2331  }
2332}
2333/* The primals are updated by the given array.
2334   Returns number of infeasibilities.
2335   After rowArray will have list of cost changes
2336*/
2337int ClpSimplexPrimal::updatePrimalsInPrimal(CoinIndexedVector *rowArray,
2338  double theta,
2339  double &objectiveChange,
2340  int valuesPass)
2341{
2342  // Cost on pivot row may change - may need to change dualIn
2343  double oldCost = 0.0;
2344  if (pivotRow_ >= 0)
2345    oldCost = cost_[sequenceOut_];
2346  //rowArray->scanAndPack();
2347  double *work = rowArray->denseVector();
2348  int number = rowArray->getNumElements();
2349  int *which = rowArray->getIndices();
2350
2351  int newNumber = 0;
2352  int pivotPosition = -1;
2353  nonLinearCost_->setChangeInCost(0.0);
2354  //printf("XX 4138 sol %g lower %g upper %g cost %g status %d\n",
2355  //   solution_[4138],lower_[4138],upper_[4138],cost_[4138],status_[4138]);
2356  // allow for case where bound+tolerance == bound
2357  //double tolerance = 0.999999*primalTolerance_;
2358  double relaxedTolerance = 1.001 * primalTolerance_;
2359  int iIndex;
2360  if (!valuesPass) {
2361    for (iIndex = 0; iIndex < number; iIndex++) {
2362
2363      int iRow = which[iIndex];
2364      double alpha = work[iIndex];
2365      work[iIndex] = 0.0;
2366      int iPivot = pivotVariable_[iRow];
2367      double change = theta * alpha;
2368      double value = solution_[iPivot] - change;
2369      solution_[iPivot] = value;
2370#ifndef NDEBUG
2371      // check if not active then okay
2372      if (!active(iRow) && (specialOptions_ & 4) == 0 && pivotRow_ != -1) {
2373        // But make sure one going out is feasible
2374        if (change > 0.0) {
2375          // going down
2376          if (value <= lower_[iPivot] + primalTolerance_) {
2377            if (iPivot == sequenceOut_ && value > lower_[iPivot] - relaxedTolerance)
2378              value = lower_[iPivot];
2379            //double difference = nonLinearCost_->setOne(iPivot, value);
2380            //assert (!difference || fabs(change) > 1.0e9);
2381          }
2382        } else {
2383          // going up
2384          if (value >= upper_[iPivot] - primalTolerance_) {
2385            if (iPivot == sequenceOut_ && value < upper_[iPivot] + relaxedTolerance)
2386              value = upper_[iPivot];
2387            //double difference = nonLinearCost_->setOne(iPivot, value);
2388            //assert (!difference || fabs(change) > 1.0e9);
2389          }
2390        }
2391      }
2392#endif
2393      if (active(iRow) || theta_ < 0.0) {
2394        clearActive(iRow);
2395        // But make sure one going out is feasible
2396        if (change > 0.0) {
2397          // going down
2398          if (value <= lower_[iPivot] + primalTolerance_) {
2399            if (iPivot == sequenceOut_ && value >= lower_[iPivot] - relaxedTolerance)
2400              value = lower_[iPivot];
2401            double difference = nonLinearCost_->setOne(iPivot, value);
2402            if (difference) {
2403              if (iRow == pivotRow_)
2404                pivotPosition = newNumber;
2405              work[newNumber] = difference;
2406              //change reduced cost on this
2407              dj_[iPivot] = -difference;
2408              which[newNumber++] = iRow;
2409            }
2410          }
2411        } else {
2412          // going up
2413          if (value >= upper_[iPivot] - primalTolerance_) {
2414            if (iPivot == sequenceOut_ && value < upper_[iPivot] + relaxedTolerance)
2415              value = upper_[iPivot];
2416            double difference = nonLinearCost_->setOne(iPivot, value);
2417            if (difference) {
2418              if (iRow == pivotRow_)
2419                pivotPosition = newNumber;
2420              work[newNumber] = difference;
2421              //change reduced cost on this
2422              dj_[iPivot] = -difference;
2423              which[newNumber++] = iRow;
2424            }
2425          }
2426        }
2427      }
2428    }
2429  } else {
2430    // values pass so look at all
2431    for (iIndex = 0; iIndex < number; iIndex++) {
2432
2433      int iRow = which[iIndex];
2434      double alpha = work[iIndex];
2435      work[iIndex] = 0.0;
2436      int iPivot = pivotVariable_[iRow];
2437      double change = theta * alpha;
2438      double value = solution_[iPivot] - change;
2439      solution_[iPivot] = value;
2440      clearActive(iRow);
2441      // But make sure one going out is feasible
2442      if (change > 0.0) {
2443        // going down
2444        if (value <= lower_[iPivot] + primalTolerance_) {
2445          if (iPivot == sequenceOut_ && value > lower_[iPivot] - relaxedTolerance)
2446            value = lower_[iPivot];
2447          double difference = nonLinearCost_->setOne(iPivot, value);
2448          if (difference) {
2449            if (iRow == pivotRow_)
2450              pivotPosition = newNumber;
2451            work[newNumber] = difference;
2452            //change reduced cost on this
2453            dj_[iPivot] = -difference;
2454            which[newNumber++] = iRow;
2455          }
2456        }
2457      } else {
2458        // going up
2459        if (value >= upper_[iPivot] - primalTolerance_) {
2460          if (iPivot == sequenceOut_ && value < upper_[iPivot] + relaxedTolerance)
2461            value = upper_[iPivot];
2462          double difference = nonLinearCost_->setOne(iPivot, value);
2463          if (difference) {
2464            if (iRow == pivotRow_)
2465              pivotPosition = newNumber;
2466            work[newNumber] = difference;
2467            //change reduced cost on this
2468            dj_[iPivot] = -difference;
2469            which[newNumber++] = iRow;
2470          }
2471        }
2472      }
2473    }
2474  }
2475  objectiveChange += nonLinearCost_->changeInCost();
2476  rowArray->setPacked();
2477#if 0
2478     rowArray->setNumElements(newNumber);
2479     rowArray->expand();
2480     if (pivotRow_ >= 0) {
2481          dualIn_ += (oldCost - cost_[sequenceOut_]);
2482          // update change vector to include pivot
2484          // and convert to packed
2485          rowArray->scanAndPack();
2486     } else {
2487          // and convert to packed
2488          rowArray->scanAndPack();
2489     }
2490#else
2491  if (pivotRow_ >= 0) {
2492    double dualIn = dualIn_ + (oldCost - cost_[sequenceOut_]);
2493    // update change vector to include pivot
2494    if (pivotPosition >= 0) {
2495      work[pivotPosition] -= dualIn;
2496    } else {
2497      work[newNumber] = -dualIn;
2498      which[newNumber++] = pivotRow_;
2499    }
2500  }
2501  rowArray->setNumElements(newNumber);
2502#endif
2503  return 0;
2504}
2505// Perturbs problem
2506void ClpSimplexPrimal::perturb(int type)
2507{
2508  if (perturbation_ > 100)
2510  if (perturbation_ == 100)
2511    perturbation_ = 50; // treat as normal
2512  int savePerturbation = perturbation_;
2513  int i;
2514  if (!numberIterations_)
2515    cleanStatus(); // make sure status okay
2516  // Make sure feasible bounds
2517  if (nonLinearCost_) {
2518    nonLinearCost_->checkInfeasibilities();
2519    //nonLinearCost_->feasibleBounds();
2520  }
2521  // look at element range
2522  double smallestNegative;
2523  double largestNegative;
2524  double smallestPositive;
2525  double largestPositive;
2526  matrix_->rangeOfElements(smallestNegative, largestNegative,
2527    smallestPositive, largestPositive);
2528  smallestPositive = CoinMin(fabs(smallestNegative), smallestPositive);
2529  largestPositive = CoinMax(fabs(largestNegative), largestPositive);
2530  double elementRatio = largestPositive / smallestPositive;
2531  if (!numberIterations_ && perturbation_ == 50) {
2532    // See if we need to perturb
2533    int numberTotal = CoinMax(numberRows_, numberColumns_);
2534    double *sort = new double[numberTotal];
2535    int nFixed = 0;
2536    for (i = 0; i < numberRows_; i++) {
2537      double lo = fabs(rowLower_[i]);
2538      double up = fabs(rowUpper_[i]);
2539      double value = 0.0;
2540      if (lo && lo < 1.0e20) {
2541        if (up && up < 1.0e20) {
2542          value = 0.5 * (lo + up);
2543          if (lo == up)
2544            nFixed++;
2545        } else {
2546          value = lo;
2547        }
2548      } else {
2549        if (up && up < 1.0e20)
2550          value = up;
2551      }
2552      sort[i] = value;
2553    }
2554    std::sort(sort, sort + numberRows_);
2555    int number = 1;
2556    double last = sort[0];
2557    for (i = 1; i < numberRows_; i++) {
2558      if (last != sort[i])
2559        number++;
2560      last = sort[i];
2561    }
2562#ifdef KEEP_GOING_IF_FIXED
2563    //printf("ratio number diff rhs %g (%d %d fixed), element ratio %g\n",((double)number)/((double) numberRows_),
2564    //   numberRows_,nFixed,elementRatio);
2565#endif
2566    if (number * 4 > numberRows_ || elementRatio > 1.0e12) {
2567      perturbation_ = 100;
2568      delete[] sort;
2569      return; // good enough
2570    }
2571    number = 0;
2572#ifdef KEEP_GOING_IF_FIXED
2573    if (!integerType_) {
2574      // look at columns
2575      nFixed = 0;
2576      for (i = 0; i < numberColumns_; i++) {
2577        double lo = fabs(columnLower_[i]);
2578        double up = fabs(columnUpper_[i]);
2579        double value = 0.0;
2580        if (lo && lo < 1.0e20) {
2581          if (up && up < 1.0e20) {
2582            value = 0.5 * (lo + up);
2583            if (lo == up)
2584              nFixed++;
2585          } else {
2586            value = lo;
2587          }
2588        } else {
2589          if (up && up < 1.0e20)
2590            value = up;
2591        }
2592        sort[i] = value;
2593      }
2594      std::sort(sort, sort + numberColumns_);
2595      number = 1;
2596      last = sort[0];
2597      for (i = 1; i < numberColumns_; i++) {
2598        if (last != sort[i])
2599          number++;
2600        last = sort[i];
2601      }
2602      //printf("cratio number diff bounds %g (%d %d fixed)\n",((double)number)/((double) numberColumns_),
2603      //     numberColumns_,nFixed);
2604    }
2605#endif
2606    delete[] sort;
2607    if (number * 4 > numberColumns_) {
2608      perturbation_ = 100;
2609      return; // good enough
2610    }
2611  }
2612  // primal perturbation
2613  double perturbation = 1.0e-20;
2614  double bias = 1.0;
2615  int numberNonZero = 0;
2616  // maximum fraction of rhs/bounds to perturb
2617  double maximumFraction = 1.0e-5;
2618  if (perturbation_ >= 50) {
2619    perturbation = 1.0e-4;
2620    for (i = 0; i < numberColumns_ + numberRows_; i++) {
2621      if (upper_[i] > lower_[i] + primalTolerance_) {
2622        double lowerValue, upperValue;
2623        if (lower_[i] > -1.0e20)
2624          lowerValue = fabs(lower_[i]);
2625        else
2626          lowerValue = 0.0;
2627        if (upper_[i] < 1.0e20)
2628          upperValue = fabs(upper_[i]);
2629        else
2630          upperValue = 0.0;
2631        double value = CoinMax(fabs(lowerValue), fabs(upperValue));
2632        value = CoinMin(value, upper_[i] - lower_[i]);
2633#if 1
2634        if (value) {
2635          perturbation += value;
2636          numberNonZero++;
2637        }
2638#else
2639        perturbation = CoinMax(perturbation, value);
2640#endif
2641      }
2642    }
2643    if (numberNonZero)
2644      perturbation /= static_cast< double >(numberNonZero);
2645    else
2646      perturbation = 1.0e-1;
2647    if (perturbation_ > 50 && perturbation_ < 55) {
2648      // reduce
2649      while (perturbation_ < 55) {
2650        perturbation_++;
2651        perturbation *= 0.25;
2652        bias *= 0.25;
2653      }
2654    } else if (perturbation_ >= 55 && perturbation_ < 60) {
2655      // increase
2656      while (perturbation_ > 55) {
2657        perturbation_--;
2658        perturbation *= 4.0;
2659      }
2660      perturbation_ = 50;
2661    }
2662  } else if (perturbation_ < 100) {
2663    perturbation = pow(10.0, perturbation_);
2664    // user is in charge
2665    maximumFraction = 1.0;
2666  }
2667  double largestZero = 0.0;
2668  double largest = 0.0;
2669  double largestPerCent = 0.0;
2670  bool printOut = (handler_->logLevel() == 63);
2671  printOut = false; //off
2672  // Check if all slack
2673  int number = 0;
2674  int iSequence;
2675  for (iSequence = 0; iSequence < numberRows_; iSequence++) {
2676    if (getRowStatus(iSequence) == basic)
2677      number++;
2678  }
2679  if (rhsScale_ > 100.0) {
2680    // tone down perturbation
2681    maximumFraction *= 0.1;
2682  }
2683  if (savePerturbation == 51) {
2684    perturbation = CoinMin(0.1, perturbation);
2685    maximumFraction *= 0.1;
2686  }
2687  if (number != numberRows_)
2688    type = 1;
2689    // modify bounds
2690    // Change so at least 1.0e-5 and no more than 0.1
2691    // For now just no more than 0.1
2692    // printf("Pert type %d perturbation %g, maxF %g\n",type,perturbation,maximumFraction);
2693    // seems much slower???#define SAVE_PERT
2694#ifdef SAVE_PERT
2695  if (2 * numberColumns_ > maximumPerturbationSize_) {
2696    delete[] perturbationArray_;
2697    maximumPerturbationSize_ = 2 * numberColumns_;
2698    perturbationArray_ = new double[maximumPerturbationSize_];
2699    for (int iColumn = 0; iColumn < maximumPerturbationSize_; iColumn++) {
2700      perturbationArray_[iColumn] = randomNumberGenerator_.randomDouble();
2701    }
2702  }
2703#endif
2704  if (type == 1) {
2705    double tolerance = 100.0 * primalTolerance_;
2706    tolerance = 10.0 * primalTolerance_; // try smaller
2707    //double multiplier = perturbation*maximumFraction;
2708    for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
2709      if (getStatus(iSequence) == basic) {
2710        double lowerValue = lower_[iSequence];
2711        double upperValue = upper_[iSequence];
2712        if (upperValue > lowerValue + tolerance) {
2713          double solutionValue = solution_[iSequence];
2714          double difference = upperValue - lowerValue;
2715          difference = CoinMin(difference, perturbation);
2716          difference = CoinMin(difference, fabs(solutionValue) + 1.0);
2717          double value = maximumFraction * (difference + bias);
2718          value = CoinMin(value, 0.1);
2719          value = CoinMax(value, primalTolerance_);
2720#ifndef SAVE_PERT
2721          value *= randomNumberGenerator_.randomDouble();
2722#else
2723          value *= perturbationArray_[2 * iSequence];
2724#endif
2725          if (value) {
2726            while (value < tolerance)
2727              value *= 3.0;
2728          }
2729          if (solutionValue - lowerValue <= primalTolerance_) {
2730            lower_[iSequence] -= value;
2731          } else if (upperValue - solutionValue <= primalTolerance_) {
2732            upper_[iSequence] += value;
2733          } else {
2734#if 0
2735                              if (iSequence >= numberColumns_) {
2736                                   // may not be at bound - but still perturb (unless free)
2737                                   if (upperValue > 1.0e30 && lowerValue < -1.0e30)
2738                                        value = 0.0;
2739                                   else
2740                                        value = - value; // as -1.0 in matrix
2741                              } else {
2742                                   value = 0.0;
2743                              }
2744#else
2745            value = 0.0;
2746#endif
2747          }
2748          if (value) {
2749            if (printOut)
2750              printf("col %d lower from %g to %g, upper from %g to %g\n",
2751                iSequence, lowerValue, lower_[iSequence], upperValue, upper_[iSequence]);
2752            if (solutionValue) {
2753              largest = CoinMax(largest, value);
2754              if (value > (fabs(solutionValue) + 1.0) * largestPerCent)
2755                largestPerCent = value / (fabs(solutionValue) + 1.0);
2756            } else {
2757              largestZero = CoinMax(largestZero, value);
2758            }
2759          }
2760        }
2761      }
2762    }
2763  } else {
2764    double tolerance = 100.0 * primalTolerance_;
2765    tolerance = 10.0 * primalTolerance_; // try smaller
2766    for (i = 0; i < numberColumns_; i++) {
2767      double lowerValue = lower_[i], upperValue = upper_[i];
2768      if (upperValue > lowerValue + primalTolerance_) {
2769        double value = perturbation * maximumFraction;
2770        value = CoinMin(value, 0.1);
2771#ifndef SAVE_PERT
2772        value *= randomNumberGenerator_.randomDouble();
2773#else
2774        value *= perturbationArray_[2 * i + 1];
2775#endif
2776        value *= randomNumberGenerator_.randomDouble();
2777        if (savePerturbation != 50) {
2778          if (fabs(value) <= primalTolerance_)
2779            value = 0.0;
2780        }
2781        if (value) {
2782          double valueL = value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
2783          // get in range
2784          if (valueL <= tolerance) {
2785            valueL *= 10.0;
2786            while (valueL <= tolerance)
2787              valueL *= 10.0;
2788          } else if (valueL > 1.0e-3) {
2789            valueL *= 0.1;
2790            while (valueL > 1.0e-3)
2791              valueL *= 0.1;
2792          }
2793          if (lowerValue > -1.0e20 && lowerValue)
2794            lowerValue -= valueL;
2795          double valueU = value * (CoinMax(1.0e-2, 1.0e-5 * fabs(upperValue)));
2796          // get in range
2797          if (valueU <= tolerance) {
2798            valueU *= 10.0;
2799            while (valueU <= tolerance)
2800              valueU *= 10.0;
2801          } else if (valueU > 1.0e-3) {
2802            valueU *= 0.1;
2803            while (valueU > 1.0e-3)
2804              valueU *= 0.1;
2805          }
2806          if (upperValue < 1.0e20 && upperValue)
2807            upperValue += valueU;
2808        }
2809        if (lowerValue != lower_[i]) {
2810          double difference = fabs(lowerValue - lower_[i]);
2811          largest = CoinMax(largest, difference);
2812          if (difference > fabs(lower_[i]) * largestPerCent)
2813            largestPerCent = fabs(difference / lower_[i]);
2814        }
2815        if (upperValue != upper_[i]) {
2816          double difference = fabs(upperValue - upper_[i]);
2817          largest = CoinMax(largest, difference);
2818          if (difference > fabs(upper_[i]) * largestPerCent)
2819            largestPerCent = fabs(difference / upper_[i]);
2820        }
2821        if (printOut)
2822          printf("col %d lower from %g to %g, upper from %g to %g\n",
2823            i, lower_[i], lowerValue, upper_[i], upperValue);
2824      }
2825      lower_[i] = lowerValue;
2826      upper_[i] = upperValue;
2827    }
2828    // biased
2829    const double *rowLower = rowLower_ - numberColumns_;
2830    const double *rowUpper = rowUpper_ - numberColumns_;
2831    for (; i < numberColumns_ + numberRows_; i++) {
2832      double lowerValue = lower_[i], upperValue = upper_[i];
2833      double value = perturbation * maximumFraction;
2834      value = CoinMin(value, 0.1);
2835      value *= randomNumberGenerator_.randomDouble();
2836      if (rowLower[i] != rowUpper[i] && upperValue > lowerValue + tolerance) {
2837        if (savePerturbation != 50) {
2838          if (fabs(value) <= primalTolerance_)
2839            value = 0.0;
2840          if (lowerValue > -1.0e20 && lowerValue)
2841            lowerValue -= value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
2842          if (upperValue < 1.0e20 && upperValue)
2843            upperValue += value * (CoinMax(1.0e-2, 1.0e-5 * fabs(upperValue)));
2844        } else if (value) {
2845          double valueL = value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
2846          // get in range
2847          if (valueL <= tolerance) {
2848            valueL *= 10.0;
2849            while (valueL <= tolerance)
2850              valueL *= 10.0;
2851          } else if (valueL > 1.0) {
2852            valueL *= 0.1;
2853            while (valueL > 1.0)
2854              valueL *= 0.1;
2855          }
2856          if (lowerValue > -1.0e20 && lowerValue)
2857            lowerValue -= valueL;
2858          double valueU = value * (CoinMax(1.0e-2, 1.0e-5 * fabs(upperValue)));
2859          // get in range
2860          if (valueU <= tolerance) {
2861            valueU *= 10.0;
2862            while (valueU <= tolerance)
2863              valueU *= 10.0;
2864          } else if (valueU > 1.0) {
2865            valueU *= 0.1;
2866            while (valueU > 1.0)
2867              valueU *= 0.1;
2868          }
2869          if (upperValue < 1.0e20 && upperValue)
2870            upperValue += valueU;
2871        }
2872      } else if (upperValue > 0.0) {
2873        //upperValue -= value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
2874        // lowerValue -= value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
2875      } else if (upperValue < 0.0) {
2876        // upperValue += value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
2877        // lowerValue += value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
2878      } else {
2879      }
2880      if (lowerValue != lower_[i]) {
2881        double difference = fabs(lowerValue - lower_[i]);
2882        largest = CoinMax(largest, difference);
2883        if (difference > fabs(lower_[i]) * largestPerCent)
2884          largestPerCent = fabs(difference / lower_[i]);
2885      }
2886      if (upperValue != upper_[i]) {
2887        double difference = fabs(upperValue - upper_[i]);
2888        largest = CoinMax(largest, difference);
2889        if (difference > fabs(upper_[i]) * largestPerCent)
2890          largestPerCent = fabs(difference / upper_[i]);
2891      }
2892      if (printOut)
2893        printf("row %d lower from %g to %g, upper from %g to %g\n",
2894          i - numberColumns_, lower_[i], lowerValue, upper_[i], upperValue);
2895      lower_[i] = lowerValue;
2896      upper_[i] = upperValue;
2897    }
2898  }
2899  // Clean up
2900  for (i = 0; i < numberColumns_ + numberRows_; i++) {
2901    switch (getStatus(i)) {
2902
2903    case basic:
2904      break;
2905    case atUpperBound:
2906      solution_[i] = upper_[i];
2907      break;
2908    case isFixed:
2909    case atLowerBound:
2910      solution_[i] = lower_[i];
2911      break;
2912    case isFree:
2913      break;
2914    case superBasic:
2915      break;
2916    }
2917  }
2918  if (largest || largestZero) {
2919    handler_->message(CLP_SIMPLEX_PERTURB, messages_)
2920      << 100.0 * maximumFraction << perturbation << largest << 100.0 * largestPerCent << largestZero
2921      << CoinMessageEol;
2922    // say perturbed
2923    perturbation_ = 101;
2924  } else {
2925    // say no need
2926    perturbation_ = 102;
2927  }
2928}
2929// un perturb
2930bool ClpSimplexPrimal::unPerturb()
2931{
2932  if (perturbation_ != 101)
2933    return false;
2934  // put back original bounds and costs
2935  createRim(1 + 4);
2936  sanityCheck();
2937  // unflag
2938  unflag();
2939  // get a valid nonlinear cost function
2940  delete nonLinearCost_;
2941  nonLinearCost_ = new ClpNonLinearCost(this);
2942  perturbation_ = 102; // stop any further perturbation
2943  // move non basic variables to new bounds
2944  nonLinearCost_->checkInfeasibilities(0.0);
2945#if 1
2946  // Try using dual
2947  return true;
2948#else
2949  gutsOfSolution(NULL, NULL, ifValuesPass != 0);
2950  return false;
2951#endif
2952}
2953// Unflag all variables and return number unflagged
2954int ClpSimplexPrimal::unflag()
2955{
2956  int i;
2957  int number = numberRows_ + numberColumns_;
2958  int numberFlagged = 0;
2959  // we can't really trust infeasibilities if there is dual error
2960  // allow tolerance bigger than standard to check on duals
2961  double relaxedToleranceD = dualTolerance_ + CoinMin(1.0e-2, 10.0 * largestDualError_);
2962  for (i = 0; i < number; i++) {
2963    if (flagged(i)) {
2964      clearFlagged(i);
2965      // only say if reasonable dj
2966      if (fabs(dj_[i]) > relaxedToleranceD)
2967        numberFlagged++;
2968    }
2969  }
2970  numberFlagged += matrix_->generalExpanded(this, 8, i);
2971  if (handler_->logLevel() > 2 && numberFlagged && objective_->type() > 1)
2972    printf("%d unflagged\n", numberFlagged);
2973  return numberFlagged;
2974}
2975// Do not change infeasibility cost and always say optimal
2976void ClpSimplexPrimal::alwaysOptimal(bool onOff)
2977{
2978  if (onOff)
2979    specialOptions_ |= 1;
2980  else
2981    specialOptions_ &= ~1;
2982}
2983bool ClpSimplexPrimal::alwaysOptimal() const
2984{
2985  return (specialOptions_ & 1) != 0;
2986}
2987// Flatten outgoing variables i.e. - always to exact bound
2988void ClpSimplexPrimal::exactOutgoing(bool onOff)
2989{
2990  if (onOff)
2991    specialOptions_ |= 4;
2992  else
2993    specialOptions_ &= ~4;
2994}
2995bool ClpSimplexPrimal::exactOutgoing() const
2996{
2997  return (specialOptions_ & 4) != 0;
2998}
2999#if ALT_UPDATE_WEIGHTS
3000static void doAlternate(ClpSimplex *model)
3001{
3002  // copy rowArray_[1]
3003  ClpPrimalColumnSteepest *steep = dynamic_cast< ClpPrimalColumnSteepest * >(model->primalColumnPivot());
3004  if (steep) {
3005#if ALT_UPDATE_WEIGHTS == 1
3006    if (!altVector[0]) {
3007      altVector[1] = new CoinIndexedVector(2000);
3008    }
3009    CoinIndexedVector *alternateWeights2 = altVector[1];
3010#else
3011    CoinIndexedVector *alternateWeights2 = steep->alternateWeights();
3012#endif
3013    double *other2 = alternateWeights2->denseVector();
3014    int *index2 = alternateWeights2->getIndices();
3015    //rowArray packed
3016    double *array = model->rowArray(1)->denseVector();
3017    int *index = model->rowArray(1)->getIndices();
3018    int n = model->rowArray(1)->getNumElements();
3019    const int *pivotVariable = model->pivotVariable();
3020    int n2 = 0;
3021    alternateWeights2->clear();
3022    if (steep->mode() != 1) {
3023      for (int i = 0; i < n; i++) {
3024        int iRow = index[i];
3025        int iPivot = pivotVariable[iRow];
3026        if (steep->reference(iPivot)) {
3027          other2[iRow] = -2.0 * array[i];
3028          index2[n2++] = iRow;
3029        }
3030      }
3031    } else {
3032      for (int i = 0; i < n; i++) {
3033        int iRow = index[i];
3034        other2[iRow] = -2.0 * array[i];
3035        index2[n2++] = iRow;
3036      }
3037    }
3038    alternateWeights2->setNumElements(n2);
3039    model->factorization()->updateColumnTranspose(model->rowArray(0),
3040      alternateWeights2);
3041  }
3042}
3043#endif
3044/*
3045  Reasons to come out (normal mode/user mode):
3046  -1 normal
3047  -2 factorize now - good iteration/ NA
3048  -3 slight inaccuracy - refactorize - iteration done/ same but factor done
3049  -4 inaccuracy - refactorize - no iteration/ NA
3050  -5 something flagged - go round again/ pivot not possible
3051  +2 looks unbounded
3052  +3 max iterations (iteration done)
3053*/
3054int ClpSimplexPrimal::pivotResult(int ifValuesPass)
3055{
3056
3057  bool roundAgain = true;
3058  int returnCode = -1;
3059
3060  // loop round if user setting and doing refactorization
3061  while (roundAgain) {
3062    roundAgain = false;
3063    returnCode = -1;
3064    pivotRow_ = -1;
3065    sequenceOut_ = -1;
3066    rowArray_[1]->clear();
3067#if 0
3068          {
3069               int seq[] = {612, 643};
3070               int k;
3071               for (k = 0; k < sizeof(seq) / sizeof(int); k++) {
3072                    int iSeq = seq[k];
3073                    if (getColumnStatus(iSeq) != basic) {
3074                         double djval;
3075                         double * work;
3076                         int number;
3077                         int * which;
3078
3079                         int iIndex;
3080                         unpack(rowArray_[1], iSeq);
3081                         factorization_->updateColumn(rowArray_[2], rowArray_[1]);
3082                         djval = cost_[iSeq];
3083                         work = rowArray_[1]->denseVector();
3084                         number = rowArray_[1]->getNumElements();
3085                         which = rowArray_[1]->getIndices();
3086
3087                         for (iIndex = 0; iIndex < number; iIndex++) {
3088
3089                              int iRow = which[iIndex];
3090                              double alpha = work[iRow];
3091                              int iPivot = pivotVariable_[iRow];
3092                              djval -= alpha * cost_[iPivot];
3093                         }
3094                         double comp = 1.0e-8 + 1.0e-7 * (CoinMax(fabs(dj_[iSeq]), fabs(djval)));
3095                         if (fabs(djval - dj_[iSeq]) > comp)
3096                              printf("Bad dj %g for %d - true is %g\n",
3097                                     dj_[iSeq], iSeq, djval);
3098                         assert (fabs(djval) < 1.0e-3 || djval * dj_[iSeq] > 0.0);
3099                         rowArray_[1]->clear();
3100                    }
3101               }
3102          }
3103#endif
3104
3105    // we found a pivot column
3106    // update the incoming column
3107    unpackPacked(rowArray_[1]);
3108    // save reduced cost
3109    double saveDj = dualIn_;
3110    factorization_->updateColumnFT(rowArray_[2], rowArray_[1]);
3111    // Get extra rows
3112    matrix_->extendUpdated(this, rowArray_[1], 0);
3113#ifdef ALT_UPDATE_WEIGHTS
3114    cilk_spawn doAlternate(this);
3115#endif
3116    // do ratio test and re-compute dj
3117#ifdef CLP_USER_DRIVEN
3118    if (solveType_ != 2 || (moreSpecialOptions_ & 512) == 0) {
3119#endif
3120      primalRow(rowArray_[1], rowArray_[3], rowArray_[2],
3121        ifValuesPass);
3122#ifdef CLP_USER_DRIVEN
3123      // user can tell which use it is
3124      int status = eventHandler_->event(ClpEventHandler::pivotRow);
3125      if (status >= 0) {
3126        problemStatus_ = 5;
3127        secondaryStatus_ = ClpEventHandler::pivotRow;
3128        break;
3129      }
3130    } else {
3131      int status = eventHandler_->event(ClpEventHandler::pivotRow);
3132      if (status >= 0) {
3133        problemStatus_ = 5;
3134        secondaryStatus_ = ClpEventHandler::pivotRow;
3135        break;
3136      }
3137    }
3138#endif
3139    if (ifValuesPass) {
3140      saveDj = dualIn_;
3141      //assert (fabs(alpha_)>=1.0e-5||(objective_->type()<2||!objective_->activated())||pivotRow_==-2);
3142      if (pivotRow_ == -1 || (pivotRow_ >= 0 && fabs(alpha_) < 1.0e-5)) {
3143        if (fabs(dualIn_) < 1.0e2 * dualTolerance_ && objective_->type() < 2) {
3144          // try other way
3145          directionIn_ = -directionIn_;
3146          primalRow(rowArray_[1], rowArray_[3], rowArray_[2],
3147            0);
3148        }
3149        if (pivotRow_ == -1 || (pivotRow_ >= 0 && fabs(alpha_) < 1.0e-5)) {
3150          if (solveType_ == 1) {
3151            // reject it
3152            char x = isColumn(sequenceIn_) ? 'C' : 'R';
3153            handler_->message(CLP_SIMPLEX_FLAG, messages_)
3154              << x << sequenceWithin(sequenceIn_)
3155              << CoinMessageEol;
3156            setFlagged(sequenceIn_);
3158            lastBadIteration_ = numberIterations_; // say be more cautious
3159            clearAll();
3160            pivotRow_ = -1;
3161          }
3162          returnCode = -5;
3163          break;
3164        }
3165      }
3166    }
3167#ifdef ALT_UPDATE_WEIGHTS
3168    cilk_sync;
3169#endif
3170    // need to clear toIndex_ in gub
3171    // ? when can I clear stuff
3172    // Clean up any gub stuff
3173    matrix_->extendUpdated(this, rowArray_[1], 1);
3174    double checkValue = 1.0e-2;
3175    if (largestDualError_ > 1.0e-5)
3176      checkValue = 1.0e-1;
3177    double test2 = dualTolerance_;
3178    double test1 = 1.0e-20;
3179#if 0 //def FEB_TRY
3180          if (factorization_->pivots() < 1) {
3181               test1 = -1.0e-4;
3182               if ((saveDj < 0.0 && dualIn_ < -1.0e-5 * dualTolerance_) ||
3183                         (saveDj > 0.0 && dualIn_ > 1.0e-5 * dualTolerance_))
3184                    test2 = 0.0; // allow through
3185          }
3186#endif
3187    if (!ifValuesPass && solveType_ == 1 && (saveDj * dualIn_ < test1 || fabs(saveDj - dualIn_) > checkValue * (1.0 + fabs(saveDj)) || fabs(dualIn_) < test2)) {
3188      if (!(saveDj * dualIn_ > 0.0 && CoinMin(fabs(saveDj), fabs(dualIn_)) > 1.0e5)) {
3189        char x = isColumn(sequenceIn_) ? 'C' : 'R';
3190        handler_->message(CLP_PRIMAL_DJ, messages_)
3191          << x << sequenceWithin(sequenceIn_)
3192          << saveDj << dualIn_
3193          << CoinMessageEol;
3194        if (lastGoodIteration_ != numberIterations_) {
3195          clearAll();
3196          pivotRow_ = -1; // say no weights update
3197          returnCode = -4;
3198          if (lastGoodIteration_ + 1 == numberIterations_) {
3199            // not looking wonderful - try cleaning bounds
3200            // put non-basics to bounds in case tolerance moved
3201            nonLinearCost_->checkInfeasibilities(0.0);
3202          }
3203          sequenceOut_ = -1;
3204          break;
3205        } else {
3206          // take on more relaxed criterion
3207          if ((saveDj * dualIn_ < test1 || fabs(saveDj - dualIn_) > 2.0e-1 * (1.0 + fabs(dualIn_)) || fabs(dualIn_) < test2) && (fabs(saveDj) > fabs(dualIn_) || saveDj * dualIn_ < 1.0e-4 || factorization_->pivots())) {
3208            // need to reject something
3209            char x = isColumn(sequenceIn_) ? 'C' : 'R';
3210            handler_->message(CLP_SIMPLEX_FLAG, messages_)
3211              << x << sequenceWithin(sequenceIn_)
3212              << CoinMessageEol;
3213            setFlagged(sequenceIn_);
3214#if 1 //def FEB_TRY
3215  // could do conditional reset of weights to get larger djs
3216            primalColumnPivot_->saveWeights(this, 6);
3217            // Make safer?
3218            double oldTolerance = factorization_->pivotTolerance();
3219            factorization_->saferTolerances(-0.99, -1.03);
3220            if (factorization_->pivotTolerance() < 1.029 * oldTolerance
3221              && oldTolerance < 0.995
3222              && !factorization_->pivots()) {
3223#ifdef CLP_USEFUL_PRINTOUT
3224              printf("Changing pivot tolerance from %g to %g and factorizing\n",
3225                oldTolerance, factorization_->pivotTolerance());
3226#endif
3227              clearAll();
3228              pivotRow_ = -1; // say no weights update
3229              returnCode = -4;
3230              if (lastGoodIteration_ + 1 == numberIterations_) {
3231                // not looking wonderful - try cleaning bounds
3232                // put non-basics to bounds in case tolerance moved
3233                nonLinearCost_->checkInfeasibilities(0.0);
3234              }
3235              sequenceOut_ = -1;
3236              break;
3237            }
3238#endif
3240            lastBadIteration_ = numberIterations_; // say be more cautious
3241            clearAll();
3242            pivotRow_ = -1;
3243            returnCode = -5;
3244            sequenceOut_ = -1;
3245            break;
3246          }
3247        }
3248      } else {
3249        //printf("%d %g %g\n",numberIterations_,saveDj,dualIn_);
3250      }
3251    }
3252    if (pivotRow_ >= 0) {
3253#ifdef CLP_USER_DRIVEN1
3254      // Got good pivot - may need to unflag stuff
3255      userChoiceWasGood(this);
3256#endif
3257      if (solveType_ >= 2 && (moreSpecialOptions_ & 512) == 0) {
3258        // **** Coding for user interface
3259        // do ray
3260        if (solveType_ == 2)
3261          primalRay(rowArray_[1]);
3262        // update duals
3263        // as packed need to find pivot row
3264        //assert (rowArray_[1]->packedMode());
3265        //int i;
3266
3267        //alpha_ = rowArray_[1]->denseVector()[pivotRow_];
3268        CoinAssert(fabs(alpha_) > 1.0e-12);
3269        double multiplier = dualIn_ / alpha_;
3270#ifndef NDEBUG
3271        rowArray_[0]->checkClear();
3272#endif
3273        rowArray_[0]->insert(pivotRow_, multiplier);
3274        factorization_->updateColumnTranspose(rowArray_[2], rowArray_[0]);
3275        // put row of tableau in rowArray[0] and columnArray[0]
3276        matrix_->transposeTimes(this, -1.0,
3277          rowArray_[0],
3278#ifdef LONG_REGION_2
3279          rowArray_[2],
3280#else
3281          columnArray_[1],
3282#endif
3283          columnArray_[0]);
3284        // update column djs
3285        int i;
3286        int *index = columnArray_[0]->getIndices();
3287        int number = columnArray_[0]->getNumElements();
3288        double *element = columnArray_[0]->denseVector();
3289        for (i = 0; i < number; i++) {
3290          int ii = index[i];
3291          dj_[ii] += element[ii];
3292          reducedCost_[ii] = dj_[ii];
3293          element[ii] = 0.0;
3294        }
3295        columnArray_[0]->setNumElements(0);
3296        // and row djs
3297        index = rowArray_[0]->getIndices();
3298        number = rowArray_[0]->getNumElements();
3299        element = rowArray_[0]->denseVector();
3300        for (i = 0; i < number; i++) {
3301          int ii = index[i];
3302          dj_[ii + numberColumns_] += element[ii];
3303          dual_[ii] = dj_[ii + numberColumns_];
3304          element[ii] = 0.0;
3305        }
3306        rowArray_[0]->setNumElements(0);
3307        // check incoming
3308        CoinAssert(fabs(dj_[sequenceIn_]) < 1.0e-1);
3309      }
3310      // if stable replace in basis
3311      // If gub or odd then alpha and pivotRow may change
3312      int updateType = 0;
3313      int updateStatus = matrix_->generalExpanded(this, 3, updateType);
3314      if (updateType >= 0)
3316          rowArray_[2],
3317          rowArray_[1],
3318          pivotRow_,
3319          alpha_,
3320          (moreSpecialOptions_ & 16) != 0);
3321
3322      // if no pivots, bad update but reasonable alpha - take and invert
3323      if (updateStatus == 2 && lastGoodIteration_ == numberIterations_ && fabs(alpha_) > 1.0e-5)
3326        // slight error
3327        if (factorization_->pivots() > 5 || updateStatus == 4) {
3328          returnCode = -3;
3329        }
3330      } else if (updateStatus == 2) {
3331        // major error
3332        // better to have small tolerance even if slower
3333        factorization_->zeroTolerance(CoinMin(factorization_->zeroTolerance(), 1.0e-15));
3334        int maxFactor = factorization_->maximumPivots();
3335        if (maxFactor > 10) {
3336          if (forceFactorization_ < 0)
3337            forceFactorization_ = maxFactor;
3338          forceFactorization_ = CoinMax(1, (forceFactorization_ >> 1));
3339        }
3340        // later we may need to unwind more e.g. fake bounds
3341        if (lastGoodIteration_ != numberIterations_) {
3342          clearAll();
3343          pivotRow_ = -1;
3344          if (solveType_ == 1 || (moreSpecialOptions_ & 512) != 0) {
3345            returnCode = -4;
3346            break;
3347          } else {
3348            // user in charge - re-factorize
3349            int lastCleaned = 0;
3350            ClpSimplexProgress dummyProgress;
3351            if (saveStatus_)
3352              statusOfProblemInPrimal(lastCleaned, 1, &dummyProgress, true, ifValuesPass);
3353            else
3354              statusOfProblemInPrimal(lastCleaned, 0, &dummyProgress, true, ifValuesPass);
3355            roundAgain = true;
3356            continue;
3357          }
3358        } else {
3359          // need to reject something
3360          if (solveType_ == 1) {
3361            char x = isColumn(sequenceIn_) ? 'C' : 'R';
3362            handler_->message(CLP_SIMPLEX_FLAG, messages_)
3363              << x << sequenceWithin(sequenceIn_)
3364              << CoinMessageEol;
3365            setFlagged(sequenceIn_);
3367          }
3368          lastBadIteration_ = numberIterations_; // say be more cautious
3369          clearAll();
3370          pivotRow_ = -1;
3371          sequenceOut_ = -1;
3372          returnCode = -5;
3373          break;
3374        }
3375      } else if (updateStatus == 3) {
3376        // out of memory
3377        // increase space if not many iterations
3378        if (factorization_->pivots() < 0.5 * factorization_->maximumPivots() && factorization_->pivots() < 200)
3379          factorization_->areaFactor(
3380            factorization_->areaFactor() * 1.1);
3381        returnCode = -2; // factorize now
3382      } else if (updateStatus == 5) {
3383        problemStatus_ = -2; // factorize now
3384      }
3385      // here do part of steepest - ready for next iteration
3386      if (!ifValuesPass)
3387        primalColumnPivot_->updateWeights(rowArray_[1]);
3388    } else {
3389      if (pivotRow_ == -1) {
3390        // no outgoing row is valid
3391        if (valueOut_ != COIN_DBL_MAX) {
3392          double objectiveChange = 0.0;
3393          theta_ = valueOut_ - valueIn_;
3394          updatePrimalsInPrimal(rowArray_[1], theta_, objectiveChange, ifValuesPass);
3395          solution_[sequenceIn_] += theta_;
3396        }
3397        rowArray_[0]->clear();
3398#ifdef CLP_USER_DRIVEN1
3399        /* Note if valueOut_ < COIN_DBL_MAX and
3400                       theta_ reasonable then this may be a valid sub flip */
3401        if (!userChoiceValid2(this)) {
3402          if (factorization_->pivots() < 5) {
3403            // flag variable
3404            char x = isColumn(sequenceIn_) ? 'C' : 'R';
3405            handler_->message(CLP_SIMPLEX_FLAG, messages_)
3406              << x << sequenceWithin(sequenceIn_)
3407              << CoinMessageEol;
3408            setFlagged(sequenceIn_);
3410            roundAgain = true;
3411            continue;
3412          } else {
3413            // try refactorizing first
3414            returnCode = 4; //say looks odd but has iterated
3415            break;
3416          }
3417        }
3418#endif
3419        if (!factorization_->pivots() && acceptablePivot_ <= 1.0e-8) {
3420          returnCode = 2; //say looks unbounded
3421          // do ray
3422          if (!nonLinearCost_->sumInfeasibilities())
3423            primalRay(rowArray_[1]);
3424        } else if (solveType_ == 2 && (moreSpecialOptions_ & 512) == 0) {
3425          // refactorize
3426          int lastCleaned = 0;
3427          ClpSimplexProgress dummyProgress;
3428          if (saveStatus_)
3429            statusOfProblemInPrimal(lastCleaned, 1, &dummyProgress, true, ifValuesPass);
3430          else
3431            statusOfProblemInPrimal(lastCleaned, 0, &dummyProgress, true, ifValuesPass);
3432          roundAgain = true;
3433          continue;
3434        } else {
3435          acceptablePivot_ = 1.0e-8;
3436          returnCode = 4; //say looks unbounded but has iterated
3437        }
3438        break;
3439      } else {
3440        // flipping from bound to bound
3441      }
3442    }
3443
3444    double oldCost = 0.0;
3445    if (sequenceOut_ >= 0)
3446      oldCost = cost_[sequenceOut_];
3447    // update primal solution
3448
3449    double objectiveChange = 0.0;
3450    // after this rowArray_[1] is not empty - used to update djs
3451    // If pivot row >= numberRows then may be gub
3452    int savePivot = pivotRow_;
3453    if (pivotRow_ >= numberRows_)
3454      pivotRow_ = -1;
3455#ifdef CLP_USER_DRIVEN
3456    if (theta_ < 0.0) {
3457      if (theta_ >= -1.0e-12)
3458        theta_ = 0.0;
3459      //else
3460      //printf("negative theta %g\n",theta_);
3461    }
3462#endif
3463    updatePrimalsInPrimal(rowArray_[1], theta_, objectiveChange, ifValuesPass);
3464    pivotRow_ = savePivot;
3465
3466    double oldValue = valueIn_;
3467    if (directionIn_ == -1) {
3468      // as if from upper bound
3469      if (sequenceIn_ != sequenceOut_) {
3470        // variable becoming basic
3471        valueIn_ -= fabs(theta_);
3472      } else {
3473        valueIn_ = lowerIn_;
3474      }
3475    } else {
3476      // as if from lower bound
3477      if (sequenceIn_ != sequenceOut_) {
3478        // variable becoming basic
3479        valueIn_ += fabs(theta_);
3480      } else {
3481        valueIn_ = upperIn_;
3482      }
3483    }
3484    objectiveChange += dualIn_ * (valueIn_ - oldValue);
3485    // outgoing
3486    if (sequenceIn_ != sequenceOut_) {
3487      if (directionOut_ > 0) {
3488        valueOut_ = lowerOut_;
3489      } else {
3490        valueOut_ = upperOut_;
3491      }
3492      if (valueOut_ < lower_[sequenceOut_] - primalTolerance_)
3493        valueOut_ = lower_[sequenceOut_] - 0.9 * primalTolerance_;
3494      else if (valueOut_ > upper_[sequenceOut_] + primalTolerance_)
3495        valueOut_ = upper_[sequenceOut_] + 0.9 * primalTolerance_;
3496      // may not be exactly at bound and bounds may have changed
3497      // Make sure outgoing looks feasible
3498      directionOut_ = nonLinearCost_->setOneOutgoing(sequenceOut_, valueOut_);
3499      // May have got inaccurate
3500      //if (oldCost!=cost_[sequenceOut_])
3501      //printf("costchange on %d from %g to %g\n",sequenceOut_,
3502      //       oldCost,cost_[sequenceOut_]);
3503      if (solveType_ < 2)
3504        dj_[sequenceOut_] = cost_[sequenceOut_] - oldCost; // normally updated next iteration
3505      solution_[sequenceOut_] = valueOut_;
3506    }
3507    // change cost and bounds on incoming if primal
3508    nonLinearCost_->setOne(sequenceIn_, valueIn_);
3509    int whatNext = housekeeping(objectiveChange);
3510    //nonLinearCost_->validate();
3511#if CLP_DEBUG > 1
3512    {
3513      double sum;
3514      int ninf = matrix_->checkFeasible(this, sum);
3515      if (ninf)
3516        printf("infeas %d\n", ninf);
3517    }
3518#endif
3519    if (whatNext == 1) {
3520      returnCode = -2; // refactorize
3521    } else if (whatNext == 2) {
3522      // maximum iterations or equivalent
3523      returnCode = 3;
3524    } else if (numberIterations_ == lastGoodIteration_ + 2 * factorization_->maximumPivots()) {
3525      // done a lot of flips - be safe
3526      returnCode = -2; // refactorize
3527    }
3528    // Check event
3529    {
3530      int status = eventHandler_->event(ClpEventHandler::endOfIteration);
3531      if (status >= 0) {
3532        problemStatus_ = 5;
3533        secondaryStatus_ = ClpEventHandler::endOfIteration;
3534        returnCode = 3;
3535      }
3536#if CLP_USER_DRIVEN
3537      int in = sequenceIn_;
3538      int out = sequenceOut_;
3539      matrix_->correctSequence(this, in, out);
3540#endif
3541    }
3542  }
3543  if ((solveType_ == 2 && (moreSpecialOptions_ & 512) == 0) && (returnCode == -2 || returnCode == -3)) {
3544    // refactorize here
3545    int lastCleaned = 0;
3546    ClpSimplexProgress dummyProgress;
3547    if (saveStatus_)
3548      statusOfProblemInPrimal(lastCleaned, 1, &dummyProgress, true, ifValuesPass);
3549    else
3550      statusOfProblemInPrimal(lastCleaned, 0, &dummyProgress, true, ifValuesPass);
3551    if (problemStatus_ == 5) {
3552      COIN_DETAIL_PRINT(printf("Singular basis\n"));
3553      problemStatus_ = -1;
3554      returnCode = 5;
3555    }
3556  }
3557#ifdef CLP_DEBUG
3558  {
3559    int i;
3560    // not [1] as may have information
3561    for (i = 0; i < 4; i++) {
3562      if (i != 1)
3563        rowArray_[i]->checkClear();
3564    }
3565    for (i = 0; i < 2; i++) {
3566      columnArray_[i]->checkClear();
3567    }
3568  }
3569#endif
3570  return returnCode;
3571}
3572// Create primal ray
3573void ClpSimplexPrimal::primalRay(CoinIndexedVector *rowArray)
3574{
3575  delete[] ray_;
3576  ray_ = new double[numberColumns_];
3577  CoinZeroN(ray_, numberColumns_);
3578  int number = rowArray->getNumElements();
3579  int *index = rowArray->getIndices();
3580  double *array = rowArray->denseVector();
3581  double way = -directionIn_;
3582  int i;
3583  double zeroTolerance = 1.0e-12;
3584  if (sequenceIn_ < numberColumns_)
3585    ray_[sequenceIn_] = directionIn_;
3586  if (!rowArray->packedMode()) {
3587    for (i = 0; i < number; i++) {
3588      int iRow = index[i];
3589      int iPivot = pivotVariable_[iRow];
3590      double arrayValue = array[iRow];
3591      if (iPivot < numberColumns_ && fabs(arrayValue) >= zeroTolerance)
3592        ray_[iPivot] = way * arrayValue;
3593    }
3594  } else {
3595    for (i = 0; i < number; i++) {
3596      int iRow = index[i];
3597      int iPivot = pivotVariable_[iRow];
3598      double arrayValue = array[i];
3599      if (iPivot < numberColumns_ && fabs(arrayValue) >= zeroTolerance)
3600        ray_[iPivot] = way * arrayValue;
3601    }
3602  }
3603}
3604/* Get next superbasic -1 if none,
3605   Normal type is 1
3606   If type is 3 then initializes sorted list
3607   if 2 uses list.
3608*/
3609int ClpSimplexPrimal::nextSuperBasic(int superBasicType,
3610  CoinIndexedVector *columnArray)
3611{
3612  int returnValue = -1;
3613  bool finished = false;
3614  while (!finished) {
3615    returnValue = firstFree_;
3616    int iColumn = firstFree_ + 1;
3617    if (superBasicType > 1) {
3618      if (superBasicType > 2) {
3619        // Initialize list
3620        // Wild guess that lower bound more natural than upper
3621        int number = 0;
3622        double *work = columnArray->denseVector();
3623        int *which = columnArray->getIndices();
3624        for (iColumn = 0; iColumn < numberRows_ + numberColumns_; iColumn++) {
3625          if (!flagged(iColumn)) {
3626            if (getStatus(iColumn) == superBasic) {
3627              if (fabs(solution_[iColumn] - lower_[iColumn]) <= primalTolerance_) {
3628                solution_[iColumn] = lower_[iColumn];
3629                setStatus(iColumn, atLowerBound);
3630              } else if (fabs(solution_[iColumn] - upper_[iColumn])
3631                <= primalTolerance_) {
3632                solution_[iColumn] = upper_[iColumn];
3633                setStatus(iColumn, atUpperBound);
3634              } else if (lower_[iColumn] < -1.0e20 && upper_[iColumn] > 1.0e20) {
3635                setStatus(iColumn, isFree);
3636                break;
3637              } else if (!flagged(iColumn)) {
3638                // put ones near bounds at end after sorting
3639                work[number] = -CoinMin(0.1 * (solution_[iColumn] - lower_[iColumn]),
3640                  upper_[iColumn] - solution_[iColumn]);
3641                which[number++] = iColumn;
3642              }
3643            }
3644          }
3645        }
3646        CoinSort_2(work, work + number, which);
3647        columnArray->setNumElements(number);
3648        CoinZeroN(work, number);
3649      }
3650      int *which = columnArray->getIndices();
3651      int number = columnArray->getNumElements();
3652      if (!number) {
3653        // finished
3654        iColumn = numberRows_ + numberColumns_;
3655        returnValue = -1;
3656      } else {
3657        number--;
3658        returnValue = which[number];
3659        iColumn = returnValue;
3660        columnArray->setNumElements(number);
3661      }
3662    } else {
3663      for (; iColumn < numberRows_ + numberColumns_; iColumn++) {
3664        if (!flagged(iColumn)) {
3665          if (getStatus(iColumn) == superBasic) {
3666            if (fabs(solution_[iColumn] - lower_[iColumn]) <= primalTolerance_) {
3667              solution_[iColumn] = lower_[iColumn];
3668              setStatus(iColumn, atLowerBound);
3669            } else if (fabs(solution_[iColumn] - upper_[iColumn])
3670              <= primalTolerance_) {
3671              solution_[iColumn] = upper_[iColumn];
3672              setStatus(iColumn, atUpperBound);
3673            } else if (lower_[iColumn] < -1.0e20 && upper_[iColumn] > 1.0e20) {
3674              setStatus(iColumn, isFree);
3675              if (fabs(dj_[iColumn]) > dualTolerance_)
3676                break;
3677            } else {
3678              break;
3679            }
3680          }
3681        }
3682      }
3683    }
3684    firstFree_ = iColumn;
3685    finished = true;
3686    if (firstFree_ == numberRows_ + numberColumns_)
3687      firstFree_ = -1;
3688    if (returnValue >= 0 && getStatus(returnValue) != superBasic && getStatus(returnValue) != isFree)
3689      finished = false; // somehow picked up odd one
3690  }
3691  return returnValue;
3692}
3693void ClpSimplexPrimal::clearAll()
3694{
3695  // Clean up any gub stuff
3696  matrix_->extendUpdated(this, rowArray_[1], 1);
3697  int number = rowArray_[1]->getNumElements();
3698  int *which = rowArray_[1]->getIndices();
3699
3700  int iIndex;
3701  for (iIndex = 0; iIndex < number; iIndex++) {
3702
3703    int iRow = which[iIndex];
3704    clearActive(iRow);
3705  }
3706  rowArray_[1]->clear();
3707  // make sure any gub sets are clean
3708  matrix_->generalExpanded(this, 11, sequenceIn_);
3709}
3710// Sort of lexicographic resolve
3711int ClpSimplexPrimal::lexSolve()
3712{
3713  algorithm_ = +1;
3714  //specialOptions_ |= 4;
3715
3716  // save data
3717  ClpDataSave data = saveData();
3718  matrix_->refresh(this); // make sure matrix okay
3719
3720  // Save so can see if doing after dual
3721  int initialStatus = problemStatus_;
3722  int initialIterations = numberIterations_;
3723  int initialNegDjs = -1;
3724  // initialize - maybe values pass and algorithm_ is +1
3725  int ifValuesPass = 0;
3726#if 0
3727     // if so - put in any superbasic costed slacks
3728     // Start can skip some things in transposeTimes
3729     specialOptions_ |= 131072;
3730     if (ifValuesPass && specialOptions_ < 0x01000000) {
3731          // Get column copy
3732          const CoinPackedMatrix * columnCopy = matrix();
3733          const int * row = columnCopy->getIndices();
3734          const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
3735          const int * columnLength = columnCopy->getVectorLengths();
3736          //const double * element = columnCopy->getElements();
3737          int n = 0;
3738          for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
3739               if (columnLength[iColumn] == 1) {
3740                    Status status = getColumnStatus(iColumn);
3741                    if (status != basic && status != isFree) {
3742                         double value = columnActivity_[iColumn];
3743                         if (fabs(value - columnLower_[iColumn]) > primalTolerance_ &&
3744                                   fabs(value - columnUpper_[iColumn]) > primalTolerance_) {
3745                              int iRow = row[columnStart[iColumn]];
3746                              if (getRowStatus(iRow) == basic) {
3747                                   setRowStatus(iRow, superBasic);
3748                                   setColumnStatus(iColumn, basic);
3749                                   n++;
3750                              }
3751                         }
3752                    }
3753               }
3754          }
3755          printf("%d costed slacks put in basis\n", n);
3756     }
3757#endif
3758  double *originalCost = NULL;
3759  double *originalLower = NULL;
3760  double *originalUpper = NULL;
3761  if (!startup(0, 0)) {
3762
3763    // Set average theta
3764    nonLinearCost_->setAverageTheta(1.0e3);
3765    int lastCleaned = 0; // last time objective or bounds cleaned up
3766
3767    // Say no pivot has occurred (for steepest edge and updates)
3768    pivotRow_ = -2;
3769
3770    // This says whether to restore things etc
3771    int factorType = 0;
3772    if (problemStatus_ < 0 && perturbation_ < 100) {
3773      perturb(0);
3774      // Can't get here if values pass
3775      assert(!ifValuesPass);
3776      gutsOfSolution(NULL, NULL);
3777      if (handler_->logLevel() > 2) {
3778        handler_->message(CLP_SIMPLEX_STATUS, messages_)
3779          << numberIterations_ << objectiveValue();
3780        handler_->printing(sumPrimalInfeasibilities_ > 0.0)
3781          << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
3782        handler_->printing(sumDualInfeasibilities_ > 0.0)
3783          << sumDualInfeasibilities_ << numberDualInfeasibilities_;
3784        handler_->printing(numberDualInfeasibilitiesWithoutFree_
3785          < numberDualInfeasibilities_)
3786          << numberDualInfeasibilitiesWithoutFree_;
3787        handler_->message() << CoinMessageEol;
3788      }
3789    }
3790    ClpSimplex *saveModel = NULL;
3791    int stopSprint = -1;
3792    int sprintPass = 0;
3793    int reasonableSprintIteration = 0;
3794    int lastSprintIteration = 0;
3795    double lastObjectiveValue = COIN_DBL_MAX;
3796    // Start check for cycles
3797    progress_.fillFromModel(this);
3798    progress_.startCheck();
3799    /*
3800            Status of problem:
3801            0 - optimal
3802            1 - infeasible
3803            2 - unbounded
3804            -1 - iterating
3805            -2 - factorization wanted
3806            -3 - redo checking without factorization
3807            -4 - looks infeasible
3808            -5 - looks unbounded
3809          */
3810    originalCost = CoinCopyOfArray(cost_, numberColumns_ + numberRows_);
3811    originalLower = CoinCopyOfArray(lower_, numberColumns_ + numberRows_);
3812    originalUpper = CoinCopyOfArray(upper_, numberColumns_ + numberRows_);
3813    while (problemStatus_ < 0) {
3814      int iRow, iColumn;
3815      // clear
3816      for (iRow = 0; iRow < 4; iRow++) {
3817        rowArray_[iRow]->clear();
3818      }
3819
3820      for (iColumn = 0; iColumn < SHORT_REGION; iColumn++) {
3821        columnArray_[iColumn]->clear();
3822      }
3823
3824      // give matrix (and model costs and bounds a chance to be
3825      // refreshed (normally null)
3826      matrix_->refresh(this);
3827      // If getting nowhere - why not give it a kick
3828#if 1
3829      if (perturbation_ < 101 && numberIterations_ > 2 * (numberRows_ + numberColumns_) && (specialOptions_ & 4) == 0
3830        && initialStatus != 10) {
3831        perturb(1);
3832        matrix_->rhsOffset(this, true, false);
3833      }
3834#endif
3835      // If we have done no iterations - special
3836      if (lastGoodIteration_ == numberIterations_ && factorType)
3837        factorType = 3;
3838      if (saveModel) {
3839        // Doing sprint
3840        if (sequenceIn_ < 0 || numberIterations_ >= stopSprint) {
3841          problemStatus_ = -1;
3842          originalModel(saveModel);
3843          saveModel = NULL;
3844          if (sequenceIn_ < 0 && numberIterations_ < reasonableSprintIteration && sprintPass > 100)
3845            primalColumnPivot_->switchOffSprint();
3846          //lastSprintIteration=numberIterations_;
3847          COIN_DETAIL_PRINT(printf("End small model\n"));
3848        }
3849      }
3850
3851      // may factorize, checks if problem finished
3852      statusOfProblemInPrimal(lastCleaned, factorType, &progress_, true, ifValuesPass, saveModel);
3853      if (initialStatus == 10) {
3854        // cleanup phase
3855        if (initialIterations != numberIterations_) {
3856          if (numberDualInfeasibilities_ > 10000 && numberDualInfeasibilities_ > 10 * initialNegDjs) {
3857            // getting worse - try perturbing
3858            if (perturbation_ < 101 && (specialOptions_ & 4) == 0) {
3859              perturb(1);
3860              matrix_->rhsOffset(this, true, false);
3861              statusOfProblemInPrimal(lastCleaned, factorType, &progress_, true, ifValuesPass, saveModel);
3862            }
3863          }
3864        } else {
3865          // save number of negative djs
3866          if (!numberPrimalInfeasibilities_)
3867            initialNegDjs = numberDualInfeasibilities_;
3868          // make sure weight won't be changed
3869          if (infeasibilityCost_ == 1.0e10)
3870            infeasibilityCost_ = 1.000001e10;
3871        }
3872      }
3873      // See if sprint says redo because of problems
3874      if (numberDualInfeasibilities_ == -776) {
3875        // Need new set of variables
3876        problemStatus_ = -1;
3877        originalModel(saveModel);
3878        saveModel = NULL;
3879        //lastSprintIteration=numberIterations_;
3880        COIN_DETAIL_PRINT(printf("End small model after\n"));
3881        statusOfProblemInPrimal(lastCleaned, factorType, &progress_, true, ifValuesPass, saveModel);
3882      }
3883      int numberSprintIterations = 0;
3884      int numberSprintColumns = primalColumnPivot_->numberSprintColumns(numberSprintIterations);
3885      if (problemStatus_ == 777) {
3886        // problems so do one pass with normal
3887        problemStatus_ = -1;
3888        originalModel(saveModel);
3889        saveModel = NULL;
3890        // Skip factorization
3891        //statusOfProblemInPrimal(lastCleaned,factorType,&progress_,false,saveModel);
3892        statusOfProblemInPrimal(lastCleaned, factorType, &progress_, true, ifValuesPass, saveModel);
3893      } else if (problemStatus_ < 0 && !saveModel && numberSprintColumns && firstFree_ < 0) {
3894        int numberSort = 0;
3895        int numberFixed = 0;
3896        int numberBasic = 0;
3897        reasonableSprintIteration = numberIterations_ + 100;
3898        int *whichColumns = new int[numberColumns_];
3899        double *weight = new double[numberColumns_];
3900        int numberNegative = 0;
3901        double sumNegative = 0.0;
3902        // now massage weight so all basic in plus good djs
3903        for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
3904          double dj = dj_[iColumn];
3905          switch (getColumnStatus(iColumn)) {
3906
3907          case basic:
3908            dj = -1.0e50;
3909            numberBasic++;
3910            break;
3911          case atUpperBound:
3912            dj = -dj;
3913            break;
3914          case isFixed:
3915            dj = 1.0e50;
3916            numberFixed++;
3917            break;
3918          case atLowerBound:
3919            dj = dj;
3920            break;
3921          case isFree:
3922            dj = -100.0 * fabs(dj);
3923            break;
3924          case superBasic:
3925            dj = -100.0 * fabs(dj);
3926            break;
3927          }
3928          if (dj < -dualTolerance_ && dj > -1.0e50) {
3929            numberNegative++;
3930            sumNegative -= dj;
3931          }
3932          weight[iColumn] = dj;
3933          whichColumns[iColumn] = iColumn;
3934        }
3935        handler_->message(CLP_SPRINT, messages_)
3936          << sprintPass << numberIterations_ - lastSprintIteration << objectiveValue() << sumNegative
3937          << numberNegative
3938          << CoinMessageEol;
3939        sprintPass++;
3940        lastSprintIteration = numberIterations_;
3941        if (objectiveValue() * optimizationDirection_ > lastObjectiveValue - 1.0e-7 && sprintPass > 5) {
3942          // switch off
3943          COIN_DETAIL_PRINT(printf("Switching off sprint\n"));
3944          primalColumnPivot_->switchOffSprint();
3945        } else {
3946          lastObjectiveValue = objectiveValue() * optimizationDirection_;
3947          // sort
3948          CoinSort_2(weight, weight + numberColumns_, whichColumns);
3949          numberSort = CoinMin(numberColumns_ - numberFixed, numberBasic + numberSprintColumns);
3950          // Sort to make consistent ?
3951          std::sort(whichColumns, whichColumns + numberSort);
3952          saveModel = new ClpSimplex(this, numberSort, whichColumns);
3953          delete[] whichColumns;
3954          delete[] weight;
3955          // Skip factorization
3956          //statusOfProblemInPrimal(lastCleaned,factorType,&progress_,false,saveModel);
3957          //statusOfProblemInPrimal(lastCleaned,factorType,&progress_,true,saveModel);
3958          stopSprint = numberIterations_ + numberSprintIterations;
3959          COIN_DETAIL_PRINT(printf("Sprint with %d columns for %d iterations\n",
3960            numberSprintColumns, numberSprintIterations));
3961        }
3962      }
3963
3964      // Say good factorization
3965      factorType = 1;
3966
3967      // Say no pivot has occurred (for steepest edge and updates)
3968      pivotRow_ = -2;
3969
3970      // exit if victory declared
3971      if (problemStatus_ >= 0) {
3972        if (originalCost) {
3973          // find number nonbasic with zero reduced costs
3974          int numberDegen = 0;
3975          int numberTotal = numberColumns_; //+numberRows_;
3976          for (int i = 0; i < numberTotal; i++) {
3977            cost_[i] = 0.0;
3978            if (getStatus(i) == atLowerBound) {
3979              if (dj_[i] <= dualTolerance_) {
3980                cost_[i] = numberTotal - i + randomNumberGenerator_.randomDouble() * 0.5;
3981                numberDegen++;
3982              } else {
3983                // fix
3984                cost_[i] = 1.0e10; //upper_[i]=lower_[i];
3985              }
3986            } else if (getStatus(i) == atUpperBound) {
3987              if (dj_[i] >= -dualTolerance_) {
3988                cost_[i] = (numberTotal - i) + randomNumberGenerator_.randomDouble() * 0.5;
3989                numberDegen++;
3990              } else {
3991                // fix
3992                cost_[i] = -1.0e10; //lower_[i]=upper_[i];
3993              }
3994            } else if (getStatus(i) == basic) {
3995              cost_[i] = (numberTotal - i) + randomNumberGenerator_.randomDouble() * 0.5;
3996            }
3997          }
3998          problemStatus_ = -1;
3999          lastObjectiveValue = COIN_DBL_MAX;
4000          // Start check for cycles
4001          progress_.fillFromModel(this);
4002          progress_.startCheck();
4003          COIN_DETAIL_PRINT(printf("%d degenerate after %d iterations\n", numberDegen,
4004            numberIterations_));
4005          if (!numberDegen) {
4006            CoinMemcpyN(originalCost, numberTotal, cost_);
4007            delete[] originalCost;
4008            originalCost = NULL;
4009            CoinMemcpyN(originalLower, numberTotal, lower_);
4010            delete[] originalLower;
4011            CoinMemcpyN(originalUpper, numberTotal, upper_);
4012            delete[] originalUpper;
4013          }
4014          delete nonLinearCost_;
4015          nonLinearCost_ = new ClpNonLinearCost(this);
4016          progress_.endOddState();
4017          continue;
4018        } else {
4019          COIN_DETAIL_PRINT(printf("exiting after %d iterations\n", numberIterations_));
4020          break;
4021        }
4022      }
4023
4024      // test for maximum iterations
4025      if (hitMaximumIterations() || (ifValuesPass == 2 && firstFree_ < 0)) {
4026        problemStatus_ = 3;
4027        break;
4028      }
4029
4030      if (firstFree_ < 0) {
4031        if (ifValuesPass) {
4032          // end of values pass
4033          ifValuesPass = 0;
4034          int status = eventHandler_->event(ClpEventHandler::endOfValuesPass);
4035          if (status >= 0) {
4036            problemStatus_ = 5;
4037            secondaryStatus_ = ClpEventHandler::endOfValuesPass;
4038            break;
4039          }
4040        }
4041      }
4042      // Check event
4043      {
4044        int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
4045        if (status >= 0) {
4046          problemStatus_ = 5;
4047          secondaryStatus_ = ClpEventHandler::endOfFactorization;
4048          break;
4049        }
4050      }
4051      // Iterate
4052      whileIterating(ifValuesPass ? 1 : 0);
4053    }
4054  }
4055  assert(!originalCost);
4056  // if infeasible get real values
4057  //printf("XXXXY final cost %g\n",infeasibilityCost_);
4058  progress_.initialWeight_ = 0.0;
4059  if (problemStatus_ == 1 && secondaryStatus_ != 6) {
4060    infeasibilityCost_ = 0.0;
4061    createRim(1 + 4);
4062    nonLinearCost_->checkInfeasibilities(0.0);
4063    sumPrimalInfeasibilities_ = nonLinearCost_->sumInfeasibilities();
4064    numberPrimalInfeasibilities_ = nonLinearCost_->numberInfeasibilities();
4065    // and get good feasible duals
4066    computeDuals(NULL);
4067  }
4068  // Stop can skip some things in transposeTimes
4069  specialOptions_ &= ~131072;
4070  // clean up
4071  unflag();
4072  finish(0);
4073  restoreData(data);
4074  return problemStatus_;
4075}
4076
4077/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
4078*/
Note: See TracBrowser for help on using the repository browser.