source: stable/1.17/Clp/src/AbcSimplexPrimal.cpp @ 2439

Last change on this file since 2439 was 2385, checked in by unxusr, 11 months ago

formatting

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