source: trunk/Clp/src/AbcSimplexDual.cpp @ 2470

Last change on this file since 2470 was 2385, checked in by unxusr, 9 months ago

formatting

  • Property svn:keywords set to Id
File size: 218.9 KB
Line 
1/* $Id: AbcSimplexDual.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 dual simplex algorithm.
7   
8   When dual feasible:
9   
10   If primal feasible, we are optimal.  Otherwise choose an infeasible
11   basic variable to leave basis (normally going to nearest bound) (B).  We
12   now need to find an incoming variable which will leave problem
13   dual feasible so we get the row of the tableau corresponding to
14   the basic variable (with the correct sign depending if basic variable
15   above or below feasibility region - as that affects whether reduced
16   cost on outgoing variable has to be positive or negative).
17   
18   We now perform a ratio test to determine which incoming variable will
19   preserve dual feasibility (C).  If no variable found then problem
20   is infeasible (in primal sense).  If there is a variable, we then
21   perform pivot and repeat.  Trivial?
22   
23   -------------------------------------------
24   
25   A) How do we get dual feasible?  If all variables have bounds then
26   it is trivial to get feasible by putting non-basic variables to
27   correct bounds.  OSL did not have a phase 1/phase 2 approach but
28   instead effectively put fake bounds on variables and this is the
29   approach here, although I had hoped to make it cleaner.
30   
31   If there is a weight of X on getting dual feasible:
32   Non-basic variables with negative reduced costs are put to
33   lesser of their upper bound and their lower bound + X.
34   Similarly, mutatis mutandis, for positive reduced costs.
35   
36   Free variables should normally be in basis, otherwise I have
37   coding which may be able to come out (and may not be correct).
38   
39   In OSL, this weight was changed heuristically, here at present
40   it is only increased if problem looks finished.  If problem is
41   feasible I check for unboundedness.  If not unbounded we
42   could play with going into primal.  As long as weights increase
43   any algorithm would be finite.
44   
45   B) Which outgoing variable to choose is a virtual base class.
46   For difficult problems steepest edge is preferred while for
47   very easy (large) problems we will need partial scan.
48   
49   C) Sounds easy, but this is hardest part of algorithm.
50   1) Instead of stopping at first choice, we may be able
51   to flip that variable to other bound and if objective
52   still improving choose again.  These mini iterations can
53   increase speed by orders of magnitude but we may need to
54   go to more of a bucket choice of variable rather than looking
55   at them one by one (for speed).
56   2) Accuracy.  Reduced costs may be of wrong sign but less than
57   tolerance.  Pivoting on these makes objective go backwards.
58   OSL modified cost so a zero move was made, Gill et al
59   (in primal analogue) modified so a strictly positive move was
60   made.  It is not quite as neat in dual but that is what we
61   try and do.  The two problems are that re-factorizations can
62   change reduced costs above and below tolerances and that when
63   finished we need to reset costs and try again.
64   3) Degeneracy.  Gill et al helps but may not be enough.  We
65   may need more.  Also it can improve speed a lot if we perturb
66   the costs significantly.
67   
68   References:
69   Forrest and Goldfarb, Steepest-edge simplex algorithms for
70   linear programming - Mathematical Programming 1992
71   Forrest and Tomlin, Implementing the simplex method for
72   the Optimization Subroutine Library - IBM Systems Journal 1992
73   Gill, Murray, Saunders, Wright A Practical Anti-Cycling
74   Procedure for Linear and Nonlinear Programming SOL report 1988
75   
76   
77   TODO:
78   
79   a) Better recovery procedures.  At present I never check on forward
80   progress.  There is checkpoint/restart with reducing
81   re-factorization frequency, but this is only on singular
82   factorizations.
83   b) Fast methods for large easy problems (and also the option for
84   the code to automatically choose which method).
85   c) We need to be able to stop in various ways for OSI - this
86   is fairly easy.
87   
88*/
89
90#include "CoinPragma.hpp"
91
92#include <math.h>
93
94#include "CoinHelperFunctions.hpp"
95#include "CoinAbcHelperFunctions.hpp"
96#include "AbcSimplexDual.hpp"
97#include "ClpEventHandler.hpp"
98#include "AbcSimplexFactorization.hpp"
99#include "AbcNonLinearCost.hpp"
100#include "CoinPackedMatrix.hpp"
101#include "CoinIndexedVector.hpp"
102#include "CoinFloatEqual.hpp"
103#include "AbcDualRowDantzig.hpp"
104#include "AbcDualRowSteepest.hpp"
105#include "ClpMessage.hpp"
106#include "ClpLinearObjective.hpp"
107class ClpSimplex;
108#include <cfloat>
109#include <cassert>
110#include <string>
111#include <stdio.h>
112#include <iostream>
113//#define CLP_DEBUG 1
114#ifdef NDEBUG
115#define NDEBUG_CLP
116#endif
117#ifndef CLP_INVESTIGATE
118#define NDEBUG_CLP
119#endif
120// dual
121
122/* *** Method
123   This is a vanilla version of dual simplex.
124   
125   It tries to be a single phase approach with a weight of 1.0 being
126   given to getting optimal and a weight of dualBound_ being
127   given to getting dual feasible.  In this version I have used the
128   idea that this weight can be thought of as a fake bound.  If the
129   distance between the lower and upper bounds on a variable is less
130   than the feasibility weight then we are always better off flipping
131   to other bound to make dual feasible.  If the distance is greater
132   then we make up a fake bound dualBound_ away from one bound.
133   If we end up optimal or primal infeasible, we check to see if
134   bounds okay.  If so we have finished, if not we increase dualBound_
135   and continue (after checking if unbounded). I am undecided about
136   free variables - there is coding but I am not sure about it.  At
137   present I put them in basis anyway.
138   
139   The code is designed to take advantage of sparsity so arrays are
140   seldom zeroed out from scratch or gone over in their entirety.
141   The only exception is a full scan to find outgoing variable.  This
142   will be changed to keep an updated list of infeasibilities (or squares
143   if steepest edge).  Also on easy problems we don't need full scan - just
144   pick first reasonable.
145   
146   One problem is how to tackle degeneracy and accuracy.  At present
147   I am using the modification of costs which I put in OSL and which was
148   extended by Gill et al.  I am still not sure of the exact details.
149   
150   The flow of dual is three while loops as follows:
151   
152   while (not finished) {
153   
154   while (not clean solution) {
155   
156   Factorize and/or clean up solution by flipping variables so
157   dual feasible.  If looks finished check fake dual bounds.
158   Repeat until status is iterating (-1) or finished (0,1,2)
159   
160   }
161   
162   while (status==-1) {
163   
164   Iterate until no pivot in or out or time to re-factorize.
165   
166   Flow is:
167   
168   choose pivot row (outgoing variable).  if none then
169   we are primal feasible so looks as if done but we need to
170   break and check bounds etc.
171   
172   Get pivot row in tableau
173   
174   Choose incoming column.  If we don't find one then we look
175   primal infeasible so break and check bounds etc.  (Also the
176   pivot tolerance is larger after any iterations so that may be
177   reason)
178   
179   If we do find incoming column, we may have to adjust costs to
180   keep going forwards (anti-degeneracy).  Check pivot will be stable
181   and if unstable throw away iteration (we will need to implement
182   flagging of basic variables sometime) and break to re-factorize.
183   If minor error re-factorize after iteration.
184   
185   Update everything (this may involve flipping variables to stay
186   dual feasible.
187   
188   }
189   
190   }
191   
192   At present we never check we are going forwards.  I overdid that in
193   OSL so will try and make a last resort.
194   
195   Needs partial scan pivot out option.
196   Needs dantzig, uninitialized and full steepest edge options (can still
197   use partial scan)
198   
199   May need other anti-degeneracy measures, especially if we try and use
200   loose tolerances as a way to solve in fewer iterations.
201   
202   I like idea of dynamic scaling.  This gives opportunity to decouple
203   different implications of scaling for accuracy, iteration count and
204   feasibility tolerance.
205   
206*/
207#define CLEAN_FIXED 0
208// Startup part of dual (may be extended to other algorithms)
209// To force to follow another run put logfile name here and define
210//#define FORCE_FOLLOW
211#ifdef FORCE_FOLLOW
212static FILE *fpFollow = NULL;
213static const char *forceFile = NULL;
214static int force_in = -1;
215static int force_out = -1;
216static int force_way_in = -1;
217static int force_way_out = -1;
218static int force_iterations = 0;
219int force_argc = 0;
220const char **force_argv = NULL;
221#endif
222void AbcSimplexDual::startupSolve()
223{
224#ifdef FORCE_FOLLOW
225  int ifld;
226  for (ifld = 1; ifld < force_argc; ifld++) {
227    if (strstr(argv[ifld], ".log")) {
228      forceFile = argv[ifld];
229      break;
230    }
231  }
232  if (!fpFollow && strlen(forceFile)) {
233    fpFollow = fopen(forceFile, "r");
234    assert(fpFollow);
235  }
236  if (fpFollow) {
237    int numberRead = atoi(argv[ifld + 1]);
238    force_iterations = atoi(argv[ifld + 1]);
239    printf("skipping %d iterations and then doing %d\n", numberRead, force_iterations);
240    for (int iteration = 0; iteration < numberRead; iteration++) {
241      // read to next Clp0102
242      char temp[300];
243      while (fgets(temp, 250, fpFollow)) {
244        if (!strncmp(temp, "dirOut", 6))
245          break;
246      }
247      sscanf(temp + 7, "%d%*s%d%*s%*s%d",
248        &force_way_out, &force_way_in);
249      fgets(temp, 250, fpFollow);
250      char cin, cout;
251      int whichIteration;
252      sscanf(temp, "%d%*f%*s%*c%c%d%*s%*c%c%d",
253        &whichIteration, &cin, &force_in, &cout, &force_out);
254      assert(whichIterations == iteration + 1);
255      if (cin == 'C')
256        force_in += numberColumns_;
257      if (cout == 'C')
258        force_out += numberColumns_;
259      printf("Iteration %d will force %d out (way %d) and %d in (way %d)\n",
260        whichIteration, force_out, force_way_out, force_in, force_way_in);
261      assert(getInternalStatus(force_out) == basic);
262      assert(getInternalStatus(force_in) != basic);
263      setInternalStatus(force_in, basic);
264      if (force_way_out == -1)
265        setInternalStatus(force_out, atUpperBound);
266      else
267        setInternalStatus(force_out, atLowerBound);
268    }
269    // clean up
270    int numberBasic = 0;
271    for (int i = 0; i < numberTotal_; i++) {
272      if (getInternalStatus(i) == basic) {
273        abcPivotVariable_[numberBasic++] = i;
274      } else if (getInternalStatus(i) == atLowerBound) {
275        if (abcLower_[i] < -1.0e30) {
276          abcLower_[i] = abcUpper_[i] - dualBound_;
277          setFakeLower(i);
278        }
279      } else if (getInternalStatus(i) == atUpperBound) {
280        if (abcUpper_[i] > 1.0e30) {
281          abcUpper_[i] = abcLower_[i] + dualBound_;
282          setFakeUpper(i);
283        }
284      }
285    }
286    assert(numberBasic == numberRows_);
287  }
288#endif
289  // initialize - no values pass and algorithm_ is -1
290  // put in standard form (and make row copy)
291  // create modifiable copies of model rim and do optional scaling
292  // If problem looks okay
293  // Do initial factorization
294  numberFake_ = 0; // Number of variables at fake bounds
295  numberChanged_ = 0; // Number of variables with changed costs
296  /// Initial coding here
297  statusOfProblemInDual(0);
298}
299void AbcSimplexDual::finishSolve()
300{
301  assert(problemStatus_ || !sumPrimalInfeasibilities_);
302}
303void AbcSimplexDual::gutsOfDual()
304{
305  double largestPrimalError = 0.0;
306  double largestDualError = 0.0;
307  lastCleaned_ = 0; // last time objective or bounds cleaned up
308  numberDisasters_ = 0;
309
310  // This says whether to restore things etc
311  //startupSolve();
312  // startup will have factorized so can skip
313  // Start check for cycles
314  abcProgress_.startCheck();
315  // Say change made on first iteration
316  changeMade_ = 1;
317  // Say last objective infinite
318  //lastObjectiveValue_=-COIN_DBL_MAX;
319  if ((stateOfProblem_ & VALUES_PASS) != 0) {
320    // modify costs so nonbasic dual feasible with fake
321    for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
322      switch (getInternalStatus(iSequence)) {
323      case basic:
324      case AbcSimplex::isFixed:
325        break;
326      case isFree:
327      case superBasic:
328        //abcCost_[iSequence]-=djSaved_[iSequence];
329        //abcDj_[iSequence]-=djSaved_[iSequence];
330        djSaved_[iSequence] = 0.0;
331        break;
332      case atLowerBound:
333        if (djSaved_[iSequence] < 0.0) {
334          //abcCost_[iSequence]-=djSaved_[iSequence];
335          //abcDj_[iSequence]-=djSaved_[iSequence];
336          djSaved_[iSequence] = 0.0;
337        }
338        break;
339      case atUpperBound:
340        if (djSaved_[iSequence] > 0.0) {
341          //abcCost_[iSequence]-=djSaved_[iSequence];
342          //abcDj_[iSequence]-=djSaved_[iSequence];
343          djSaved_[iSequence] = 0.0;
344        }
345        break;
346      }
347    }
348  }
349  progressFlag_ = 0;
350  /*
351    Status of problem:
352    0 - optimal
353    1 - infeasible
354    2 - unbounded
355    -1 - iterating
356    -2 - factorization wanted
357    -3 - redo checking without factorization
358    -4 - looks infeasible
359  */
360  int factorizationType = 1;
361  double pivotTolerance = abcFactorization_->pivotTolerance();
362  while (problemStatus_ < 0) {
363    // If getting nowhere - return
364    double limit = numberRows_ + numberColumns_;
365    limit = 10000.0 + 2000.0 * limit;
366#if ABC_NORMAL_DEBUG > 0
367    if (numberIterations_ > limit) {
368      printf("Set flag so ClpSimplexDual done\n");
369      abort();
370    }
371#endif
372    bool disaster = false;
373    if (disasterArea_ && disasterArea_->check()) {
374      disasterArea_->saveInfo();
375      disaster = true;
376    }
377    if (numberIterations_ > baseIteration_) {
378      // may factorize, checks if problem finished
379      statusOfProblemInDual(factorizationType);
380      factorizationType = 1;
381      largestPrimalError = CoinMax(largestPrimalError, largestPrimalError_);
382      largestDualError = CoinMax(largestDualError, largestDualError_);
383      if (disaster)
384        problemStatus_ = 3;
385    }
386    int saveNumberIterations = numberIterations_;
387    // exit if victory declared
388    if (problemStatus_ >= 0)
389      break;
390
391    // test for maximum iterations
392    if (hitMaximumIterations()) {
393      problemStatus_ = 3;
394      break;
395    }
396    // Check event
397    int eventStatus = eventHandler_->event(ClpEventHandler::endOfFactorization);
398    if (eventStatus >= 0) {
399      problemStatus_ = 5;
400      secondaryStatus_ = ClpEventHandler::endOfFactorization;
401      break;
402    }
403    if (!numberPrimalInfeasibilities_ && problemStatus_ < 0) {
404#if ABC_NORMAL_DEBUG > 0
405      printf("Primal feasible - sum dual %g - go to primal\n", sumDualInfeasibilities_);
406#endif
407      problemStatus_ = 10;
408      break;
409    }
410    // Do iterations ? loop round
411#if PARTITION_ROW_COPY == 1
412    stateOfProblem_ |= NEED_BASIS_SORT;
413#endif
414    while (true) {
415#if PARTITION_ROW_COPY == 1
416      if ((stateOfProblem_ & NEED_BASIS_SORT) != 0) {
417        int iVector = getAvailableArray();
418        abcMatrix_->sortUseful(usefulArray_[iVector]);
419        setAvailableArray(iVector);
420        stateOfProblem_ &= ~NEED_BASIS_SORT;
421      }
422#endif
423      problemStatus_ = -1;
424      if ((stateOfProblem_ & VALUES_PASS) != 0) {
425        assert(djSaved_ - abcDj_ == numberTotal_);
426        abcDj_ = djSaved_;
427        djSaved_ = NULL;
428      }
429#if ABC_PARALLEL
430      if (parallelMode_ == 0) {
431#endif
432        whileIteratingSerial();
433#if ABC_PARALLEL
434      } else {
435#if ABC_PARALLEL == 1
436        whileIteratingThread();
437#else
438        whileIteratingCilk();
439#endif
440      }
441#endif
442      if ((stateOfProblem_ & VALUES_PASS) != 0) {
443        djSaved_ = abcDj_;
444        abcDj_ = djSaved_ - numberTotal_;
445      }
446      if (pivotRow_ != -100) {
447        if (numberIterations_ > saveNumberIterations || problemStatus_ >= 0) {
448          if (problemStatus_ == 10 && abcFactorization_->pivots() && numberTimesOptimal_ < 3) {
449            factorizationType = 5;
450            problemStatus_ = -4;
451            if (!numberPrimalInfeasibilities_)
452              factorizationType = 9;
453          } else if (problemStatus_ == -2) {
454            factorizationType = 5;
455          }
456          break;
457        } else if (abcFactorization_->pivots() || problemStatus_ == -6 || pivotTolerance < abcFactorization_->pivotTolerance()) {
458          factorizationType = 5;
459          if (pivotTolerance < abcFactorization_->pivotTolerance()) {
460            pivotTolerance = abcFactorization_->pivotTolerance();
461            if (!abcFactorization_->pivots())
462              abcFactorization_->minimumPivotTolerance(CoinMin(0.25, 1.05 * abcFactorization_->minimumPivotTolerance()));
463          }
464          break;
465        } else {
466          bool tryFake = numberAtFakeBound() > 0 && currentDualBound_ < 1.0e17;
467          if (problemStatus_ == -5 || tryFake) {
468            if (numberTimesOptimal_ > 10)
469              problemStatus_ = 4;
470            numberTimesOptimal_++;
471            int numberChanged = -1;
472            if (numberFlagged_) {
473              numberFlagged_ = 0;
474              for (int iRow = 0; iRow < numberRows_; iRow++) {
475                int iPivot = abcPivotVariable_[iRow];
476                if (flagged(iPivot)) {
477                  clearFlagged(iPivot);
478                }
479              }
480              numberChanged = 0;
481            } else if (tryFake) {
482              double dummyChange;
483              numberChanged = changeBounds(0, dummyChange);
484              if (numberChanged > 0) {
485                handler_->message(CLP_DUAL_ORIGINAL, messages_)
486                  << CoinMessageEol;
487                bounceTolerances(3 /*(numberChanged>0) ? 3 : 1*/);
488                // temp
489#if ABC_NORMAL_DEBUG > 0
490                //printf("should not need to break out of loop\n");
491#endif
492                abcProgress_.reset();
493                break;
494              } else if (numberChanged < -1) {
495                // some flipped
496                bounceTolerances(3);
497                abcProgress_.reset();
498              }
499            } else if (numberChanged < 0) {
500              // what now
501              printf("No iterations since last inverts A - nothing done - think\n");
502              abort();
503            }
504          } else if (problemStatus_ != -4) {
505            // what now
506            printf("No iterations since last inverts B - nothing done - think\n");
507            abort();
508          }
509        }
510      } else {
511        // end of values pass
512        // we need to tell statusOf.. to restore costs etc
513        // and take off flag
514        bounceTolerances(2);
515        assert(!solution_);
516        delete[] solution_;
517        solution_ = NULL;
518        gutsOfSolution(3);
519        makeNonFreeVariablesDualFeasible();
520        stateOfProblem_ &= ~VALUES_PASS;
521        problemStatus_ = -2;
522        factorizationType = 5;
523        break;
524      }
525    }
526    if (problemStatus_ == 1 && (progressFlag_ & 8) != 0 && fabs(rawObjectiveValue_) > 1.0e10) {
527#if ABC_NORMAL_DEBUG > 0
528      printf("fix looks infeasible when has been\n"); // goto primal
529#endif
530      problemStatus_ = 10; // infeasible - but has looked feasible
531    }
532    if (!problemStatus_ && abcFactorization_->pivots()) {
533      gutsOfSolution(1); // need to compute duals
534      //computeInternalObjectiveValue();
535      checkCutoff(true);
536    }
537  }
538  largestPrimalError_ = largestPrimalError;
539  largestDualError_ = largestDualError;
540}
541/// The duals are updated by the given arrays.
542static void updateDualsInDualBit2(CoinPartitionedVector &row,
543  double *COIN_RESTRICT djs,
544  double theta, int iBlock)
545{
546  int offset = row.startPartition(iBlock);
547  double *COIN_RESTRICT work = row.denseVector() + offset;
548  const int *COIN_RESTRICT which = row.getIndices() + offset;
549  int number = row.getNumElements(iBlock);
550  for (int i = 0; i < number; i++) {
551    int iSequence = which[i];
552    double alphaI = work[i];
553    work[i] = 0.0;
554    double value = djs[iSequence] - theta * alphaI;
555    djs[iSequence] = value;
556  }
557  row.setNumElementsPartition(iBlock, 0);
558}
559static void updateDualsInDualBit(CoinPartitionedVector &row,
560  double *djs,
561  double theta, int numberBlocks)
562{
563  for (int iBlock = 1; iBlock < numberBlocks; iBlock++) {
564    cilk_spawn updateDualsInDualBit2(row, djs, theta, iBlock);
565  }
566  updateDualsInDualBit2(row, djs, theta, 0);
567  cilk_sync;
568}
569/// The duals are updated by the given arrays.
570
571void AbcSimplexDual::updateDualsInDual()
572{
573  //dual->usefulArray(3)->compact();
574  CoinPartitionedVector &array = usefulArray_[arrayForTableauRow_];
575  //array.compact();
576  int numberBlocks = array.getNumPartitions();
577  int number = array.getNumElements();
578  if (numberBlocks) {
579    if (number > 100) {
580      updateDualsInDualBit(array, abcDj_, theta_, numberBlocks);
581    } else {
582      for (int iBlock = 0; iBlock < numberBlocks; iBlock++) {
583        updateDualsInDualBit2(array, abcDj_, theta_, iBlock);
584      }
585    }
586    array.compact();
587  } else {
588    double *COIN_RESTRICT work = array.denseVector();
589    const int *COIN_RESTRICT which = array.getIndices();
590    double *COIN_RESTRICT reducedCost = abcDj_;
591#pragma cilk grainsize = 128
592    cilk_for(int i = 0; i < number; i++)
593    {
594      //for (int i = 0; i < number; i++) {
595      int iSequence = which[i];
596      double alphaI = work[i];
597      work[i] = 0.0;
598      double value = reducedCost[iSequence] - theta_ * alphaI;
599      reducedCost[iSequence] = value;
600    }
601    array.setNumElements(0);
602  }
603}
604int AbcSimplexDual::nextSuperBasic()
605{
606  if (firstFree_ >= 0) {
607    // make sure valid
608    int returnValue = firstFree_;
609    int iColumn = firstFree_ + 1;
610#ifdef PAN
611    if (fakeSuperBasic(firstFree_) < 0) {
612      // somehow cleaned up
613      returnValue = -1;
614      for (; iColumn < numberTotal_; iColumn++) {
615        if (getInternalStatus(iColumn) == isFree || getInternalStatus(iColumn) == superBasic) {
616          if (fakeSuperBasic(iColumn) >= 0) {
617            if (fabs(abcDj_[iColumn]) > 1.0e2 * dualTolerance_)
618              break;
619          }
620        }
621      }
622      if (iColumn < numberTotal_) {
623        returnValue = iColumn;
624        iColumn++;
625      }
626    }
627#endif
628    for (; iColumn < numberTotal_; iColumn++) {
629      if (getInternalStatus(iColumn) == isFree || getInternalStatus(iColumn) == superBasic) {
630#ifdef PAN
631        if (fakeSuperBasic(iColumn) >= 0) {
632#endif
633          if (fabs(abcDj_[iColumn]) > dualTolerance_)
634            break;
635#ifdef PAN
636        }
637#endif
638      }
639    }
640    firstFree_ = iColumn;
641    if (firstFree_ == numberTotal_)
642      firstFree_ = -1;
643    return returnValue;
644  } else {
645    return -1;
646  }
647}
648/*
649  Chooses dual pivot row
650  For easy problems we can just choose one of the first rows we look at
651*/
652void AbcSimplexDual::dualPivotRow()
653{
654  //int lastPivotRow=pivotRow_;
655  pivotRow_ = -1;
656  const double *COIN_RESTRICT lowerBasic = lowerBasic_;
657  const double *COIN_RESTRICT upperBasic = upperBasic_;
658  const double *COIN_RESTRICT solutionBasic = solutionBasic_;
659  const int *COIN_RESTRICT pivotVariable = abcPivotVariable_;
660  freeSequenceIn_ = -1;
661  double worstValuesPass = 0.0;
662  if ((stateOfProblem_ & VALUES_PASS) != 0) {
663    for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
664      switch (getInternalStatus(iSequence)) {
665      case basic:
666#if ABC_NORMAL_DEBUG > 0
667        if (fabs(abcDj_[iSequence]) > 1.0e-6)
668          printf("valuespass basic dj %d status %d dj %g\n", iSequence,
669            getInternalStatus(iSequence), abcDj_[iSequence]);
670#endif
671        break;
672      case AbcSimplex::isFixed:
673        break;
674      case isFree:
675      case superBasic:
676#if ABC_NORMAL_DEBUG > 0
677        if (fabs(abcDj_[iSequence]) > 1.0e-6)
678          printf("Bad valuespass dj %d status %d dj %g\n", iSequence,
679            getInternalStatus(iSequence), abcDj_[iSequence]);
680#endif
681        break;
682      case atLowerBound:
683#if ABC_NORMAL_DEBUG > 0
684        if (abcDj_[iSequence] < -1.0e-6)
685          printf("Bad valuespass dj %d status %d dj %g\n", iSequence,
686            getInternalStatus(iSequence), abcDj_[iSequence]);
687#endif
688        break;
689      case atUpperBound:
690#if ABC_NORMAL_DEBUG > 0
691        if (abcDj_[iSequence] > 1.0e-6)
692          printf("Bad valuespass dj %d status %d dj %g\n", iSequence,
693            getInternalStatus(iSequence), abcDj_[iSequence]);
694#endif
695        break;
696      }
697    }
698    // get worst basic dual - later do faster using steepest infeasibility array
699    // does not have to be worst - just bad
700    int iWorst = -1;
701    double worst = dualTolerance_;
702    for (int i = 0; i < numberRows_; i++) {
703      int iSequence = abcPivotVariable_[i];
704      if (fabs(abcDj_[iSequence]) > worst) {
705        iWorst = i;
706        worst = fabs(abcDj_[iSequence]);
707      }
708    }
709    if (iWorst >= 0) {
710      pivotRow_ = iWorst;
711      worstValuesPass = abcDj_[abcPivotVariable_[iWorst]];
712    } else {
713      // factorize and clean up costs
714      pivotRow_ = -100;
715      return;
716    }
717  }
718  if (false && numberFreeNonBasic_ > abcFactorization_->maximumPivots()) {
719    lastFirstFree_ = -1;
720    firstFree_ = -1;
721  }
722  if (firstFree_ >= 0) {
723    // free
724    int nextFree = nextSuperBasic();
725    if (nextFree >= 0) {
726      freeSequenceIn_ = nextFree;
727      // **** should skip price (and maybe weights update)
728      // unpack vector and find a good pivot
729      unpack(usefulArray_[arrayForBtran_], nextFree);
730      abcFactorization_->updateColumn(usefulArray_[arrayForBtran_]);
731
732      double *COIN_RESTRICT work = usefulArray_[arrayForBtran_].denseVector();
733      int number = usefulArray_[arrayForBtran_].getNumElements();
734      int *COIN_RESTRICT which = usefulArray_[arrayForBtran_].getIndices();
735      double bestFeasibleAlpha = 0.0;
736      int bestFeasibleRow = -1;
737      double bestInfeasibleAlpha = 0.0;
738      int bestInfeasibleRow = -1;
739      // Go for big infeasibilities unless big changes
740      double maxInfeas = (fabs(dualOut_) > 1.0e3) ? 0.001 : 1.0e6;
741      for (int i = 0; i < number; i++) {
742        int iRow = which[i];
743        double alpha = fabs(work[iRow]);
744        if (alpha > 1.0e-3) {
745          double value = solutionBasic[iRow];
746          double lower = lowerBasic[iRow];
747          double upper = upperBasic[iRow];
748          double infeasibility = 0.0;
749          if (value > upper)
750            infeasibility = value - upper;
751          else if (value < lower)
752            infeasibility = lower - value;
753          infeasibility = CoinMin(maxInfeas, infeasibility);
754          if (infeasibility * alpha > bestInfeasibleAlpha && alpha > 1.0e-1) {
755            if (!flagged(pivotVariable[iRow])) {
756              bestInfeasibleAlpha = infeasibility * alpha;
757              bestInfeasibleRow = iRow;
758            }
759          }
760          if (alpha > bestFeasibleAlpha && (lower > -1.0e20 || upper < 1.0e20)) {
761            bestFeasibleAlpha = alpha;
762            bestFeasibleRow = iRow;
763          }
764        }
765      }
766      if (bestInfeasibleRow >= 0)
767        pivotRow_ = bestInfeasibleRow;
768      else if (bestFeasibleAlpha > 1.0e-2)
769        pivotRow_ = bestFeasibleRow;
770      usefulArray_[arrayForBtran_].clear();
771      if (pivotRow_ < 0) {
772        // no good
773        freeSequenceIn_ = -1;
774        if (!abcFactorization_->pivots()) {
775          lastFirstFree_ = -1;
776          firstFree_ = -1;
777        }
778      }
779    }
780  }
781  if (pivotRow_ < 0) {
782    // put back so can test
783    //pivotRow_=lastPivotRow;
784    // get pivot row using whichever method it is
785    pivotRow_ = abcDualRowPivot_->pivotRow();
786  }
787
788  if (pivotRow_ >= 0) {
789    sequenceOut_ = pivotVariable[pivotRow_];
790    valueOut_ = solutionBasic[pivotRow_];
791    lowerOut_ = lowerBasic[pivotRow_];
792    upperOut_ = upperBasic[pivotRow_];
793    // if we have problems we could try other way and hope we get a
794    // zero pivot?
795    if (valueOut_ > upperOut_) {
796      directionOut_ = -1;
797      dualOut_ = valueOut_ - upperOut_;
798    } else if (valueOut_ < lowerOut_) {
799      directionOut_ = 1;
800      dualOut_ = lowerOut_ - valueOut_;
801    } else {
802      // odd (could be free) - it's feasible - go to nearest
803      if (valueOut_ - lowerOut_ < upperOut_ - valueOut_) {
804        directionOut_ = 1;
805        dualOut_ = lowerOut_ - valueOut_;
806      } else {
807        directionOut_ = -1;
808        dualOut_ = valueOut_ - upperOut_;
809      }
810    }
811    if (worstValuesPass) {
812      // in values pass so just use sign of dj
813      // We don't want to go through any barriers so set dualOut low
814      // free variables will never be here
815      // need to think about tolerances with inexact solutions
816      // could try crashing in variables away from bounds and with near zero djs
817      // pivot to take out bad slack and increase objective
818      double fakeValueOut = solution_[sequenceOut_];
819      dualOut_ = 1.0e-10;
820      if (upperOut_ > lowerOut_ || fabs(fakeValueOut - upperOut_) < 100.0 * primalTolerance_) {
821        if (abcDj_[sequenceOut_] > 0.0) {
822          directionOut_ = 1;
823        } else {
824          directionOut_ = -1;
825        }
826      } else {
827        // foolishly trust primal solution
828        if (fakeValueOut > upperOut_) {
829          directionOut_ = -1;
830        } else {
831          directionOut_ = 1;
832        }
833      }
834#if 1
835      // make sure plausible but only when fake primals stored
836      if (directionOut_ < 0 && fabs(fakeValueOut - upperOut_) > dualBound_ + primalTolerance_) {
837        if (fabs(fakeValueOut - upperOut_) > fabs(fakeValueOut - lowerOut_))
838          directionOut_ = 1;
839      } else if (directionOut_ > 0 && fabs(fakeValueOut - lowerOut_) > dualBound_ + primalTolerance_) {
840        if (fabs(fakeValueOut - upperOut_) < fabs(fakeValueOut - lowerOut_))
841          directionOut_ = -1;
842      }
843#endif
844    }
845#if ABC_NORMAL_DEBUG > 3
846    assert(dualOut_ >= 0.0);
847#endif
848  } else {
849    sequenceOut_ = -1;
850  }
851  alpha_ = 0.0;
852  upperTheta_ = COIN_DBL_MAX;
853  return;
854}
855// Checks if any fake bounds active - if so returns number and modifies
856// dualBound_ and everything.
857// Free variables will be left as free
858// Returns number of bounds changed if >=0
859// Returns -1 if not initialize and no effect
860// Fills in changeVector which can be used to see if unbounded
861// initialize is 0 or 3
862// and cost of change vector
863int AbcSimplexDual::changeBounds(int initialize,
864  double &changeCost)
865{
866  assert(initialize == 0 || initialize == 1 || initialize == 3 || initialize == 4);
867  numberFake_ = 0;
868  double *COIN_RESTRICT abcLower = abcLower_;
869  double *COIN_RESTRICT abcUpper = abcUpper_;
870  double *COIN_RESTRICT abcSolution = abcSolution_;
871  double *COIN_RESTRICT abcCost = abcCost_;
872  if (!initialize) {
873    int numberInfeasibilities;
874    double newBound;
875    newBound = 5.0 * currentDualBound_;
876    numberInfeasibilities = 0;
877    changeCost = 0.0;
878    // put back original bounds and then check
879    CoinAbcMemcpy(abcLower, lowerSaved_, numberTotal_);
880    CoinAbcMemcpy(abcUpper, upperSaved_, numberTotal_);
881    const double *COIN_RESTRICT dj = abcDj_;
882    int iSequence;
883    // bounds will get bigger - just look at ones at bounds
884    int numberFlipped = 0;
885    double checkBound = 1.000000000001 * currentDualBound_;
886    for (iSequence = 0; iSequence < numberTotal_; iSequence++) {
887      double lowerValue = abcLower[iSequence];
888      double upperValue = abcUpper[iSequence];
889      double value = abcSolution[iSequence];
890      setFakeBound(iSequence, AbcSimplexDual::noFake);
891      switch (getInternalStatus(iSequence)) {
892
893      case basic:
894      case AbcSimplex::isFixed:
895        break;
896      case isFree:
897      case superBasic:
898        break;
899      case atUpperBound:
900        if (fabs(value - upperValue) > primalTolerance_) {
901          // can we flip
902          if (value - lowerValue < checkBound && dj[iSequence] > -currentDualTolerance_) {
903            abcSolution_[iSequence] = lowerValue;
904            setInternalStatus(iSequence, atLowerBound);
905            numberFlipped++;
906          } else {
907            numberInfeasibilities++;
908          }
909        }
910        break;
911      case atLowerBound:
912        if (fabs(value - lowerValue) > primalTolerance_) {
913          // can we flip
914          if (upperValue - value < checkBound && dj[iSequence] < currentDualTolerance_) {
915            abcSolution_[iSequence] = upperValue;
916            setInternalStatus(iSequence, atUpperBound);
917            numberFlipped++;
918          } else {
919            numberInfeasibilities++;
920          }
921        }
922        break;
923      }
924    }
925    // If dual infeasible then carry on
926    if (numberInfeasibilities) {
927      handler_->message(CLP_DUAL_CHECKB, messages_)
928        << newBound
929        << CoinMessageEol;
930      int iSequence;
931      for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
932        double lowerValue = abcLower[iSequence];
933        double upperValue = abcUpper[iSequence];
934        double newLowerValue;
935        double newUpperValue;
936        Status status = getInternalStatus(iSequence);
937        if (status == atUpperBound || status == atLowerBound) {
938          double value = abcSolution[iSequence];
939          if (value - lowerValue <= upperValue - value) {
940            newLowerValue = CoinMax(lowerValue, value - 0.666667 * newBound);
941            newUpperValue = CoinMin(upperValue, newLowerValue + newBound);
942          } else {
943            newUpperValue = CoinMin(upperValue, value + 0.666667 * newBound);
944            newLowerValue = CoinMax(lowerValue, newUpperValue - newBound);
945          }
946          abcLower[iSequence] = newLowerValue;
947          abcUpper[iSequence] = newUpperValue;
948          if (newLowerValue > lowerValue) {
949            if (newUpperValue < upperValue) {
950              setFakeBound(iSequence, AbcSimplexDual::bothFake);
951#ifdef CLP_INVESTIGATE
952              abort(); // No idea what should happen here - I have never got here
953#endif
954              numberFake_++;
955            } else {
956              setFakeBound(iSequence, AbcSimplexDual::lowerFake);
957              assert(abcLower[iSequence] > lowerSaved_[iSequence]);
958              assert(abcUpper[iSequence] == upperSaved_[iSequence]);
959              numberFake_++;
960            }
961          } else {
962            if (newUpperValue < upperValue) {
963              setFakeBound(iSequence, AbcSimplexDual::upperFake);
964              assert(abcLower[iSequence] == lowerSaved_[iSequence]);
965              assert(abcUpper[iSequence] < upperSaved_[iSequence]);
966              numberFake_++;
967            }
968          }
969          if (status == atUpperBound)
970            abcSolution[iSequence] = newUpperValue;
971          else
972            abcSolution[iSequence] = newLowerValue;
973          double movement = abcSolution[iSequence] - value;
974          changeCost += movement * abcCost[iSequence];
975        }
976      }
977      currentDualBound_ = newBound;
978#if ABC_NORMAL_DEBUG > 0
979      printf("new dual bound %g\n", currentDualBound_);
980#endif
981    } else {
982      numberInfeasibilities = -1 - numberFlipped;
983    }
984    return numberInfeasibilities;
985  } else {
986    int iSequence;
987    if (initialize == 3) {
988      for (iSequence = 0; iSequence < numberTotal_; iSequence++) {
989#if ABC_NORMAL_DEBUG > 1
990        int bad = -1;
991        if (getFakeBound(iSequence) == noFake) {
992          if (abcLower_[iSequence] != lowerSaved_[iSequence] || abcUpper_[iSequence] != upperSaved_[iSequence])
993            bad = 0;
994        } else if (getFakeBound(iSequence) == lowerFake) {
995          if (abcLower_[iSequence] == lowerSaved_[iSequence] || abcUpper_[iSequence] != upperSaved_[iSequence])
996            bad = 1;
997        } else if (getFakeBound(iSequence) == upperFake) {
998          if (abcLower_[iSequence] != lowerSaved_[iSequence] || abcUpper_[iSequence] == upperSaved_[iSequence])
999            bad = 2;
1000        } else {
1001          if (abcLower_[iSequence] == lowerSaved_[iSequence] || abcUpper_[iSequence] == upperSaved_[iSequence])
1002            bad = 3;
1003        }
1004        if (bad >= 0)
1005          printf("%d status %d l %g,%g u %g,%g\n", iSequence, internalStatus_[iSequence],
1006            abcLower_[iSequence], lowerSaved_[iSequence],
1007            abcUpper_[iSequence], upperSaved_[iSequence]);
1008#endif
1009        setFakeBound(iSequence, AbcSimplexDual::noFake);
1010      }
1011      CoinAbcMemcpy(abcLower_, lowerSaved_, numberTotal_);
1012      CoinAbcMemcpy(abcUpper_, upperSaved_, numberTotal_);
1013    }
1014    double testBound = /*0.999999 **/ currentDualBound_;
1015    for (iSequence = 0; iSequence < numberTotal_; iSequence++) {
1016      Status status = getInternalStatus(iSequence);
1017      if (status == atUpperBound || status == atLowerBound) {
1018        double lowerValue = abcLower[iSequence];
1019        double upperValue = abcUpper[iSequence];
1020        double value = abcSolution[iSequence];
1021        if (lowerValue > -largeValue_ || upperValue < largeValue_) {
1022          if (true || lowerValue - value > -0.5 * currentDualBound_ || upperValue - value < 0.5 * currentDualBound_) {
1023            if (fabs(lowerValue - value) <= fabs(upperValue - value)) {
1024              if (upperValue > lowerValue + testBound) {
1025                if (getFakeBound(iSequence) == AbcSimplexDual::noFake)
1026                  numberFake_++;
1027                abcUpper[iSequence] = lowerValue + currentDualBound_;
1028                double difference = upperSaved_[iSequence] - abcUpper[iSequence];
1029                if (difference < 1.0e-10 * (1.0 + fabs(upperSaved_[iSequence]) + dualBound_))
1030                  abcUpper[iSequence] = upperSaved_[iSequence];
1031                assert(abcUpper[iSequence] <= upperSaved_[iSequence] + 1.0e-3 * currentDualBound_);
1032                assert(abcLower[iSequence] >= lowerSaved_[iSequence] - 1.0e-3 * currentDualBound_);
1033                if (upperSaved_[iSequence] - abcUpper[iSequence] > abcLower[iSequence] - lowerSaved_[iSequence]) {
1034                  setFakeBound(iSequence, AbcSimplexDual::upperFake);
1035                  assert(abcLower[iSequence] == lowerSaved_[iSequence]);
1036                  assert(abcUpper[iSequence] < upperSaved_[iSequence]);
1037                } else {
1038                  setFakeBound(iSequence, AbcSimplexDual::lowerFake);
1039                  assert(abcLower[iSequence] > lowerSaved_[iSequence]);
1040                  assert(abcUpper[iSequence] == upperSaved_[iSequence]);
1041                }
1042              }
1043            } else {
1044              if (lowerValue < upperValue - testBound) {
1045                if (getFakeBound(iSequence) == AbcSimplexDual::noFake)
1046                  numberFake_++;
1047                abcLower[iSequence] = upperValue - currentDualBound_;
1048                double difference = abcLower[iSequence] - lowerSaved_[iSequence];
1049                if (difference < 1.0e-10 * (1.0 + fabs(lowerSaved_[iSequence]) + dualBound_))
1050                  abcLower[iSequence] = lowerSaved_[iSequence];
1051                assert(abcUpper[iSequence] <= upperSaved_[iSequence] + 1.0e-3 * currentDualBound_);
1052                assert(abcLower[iSequence] >= lowerSaved_[iSequence] - 1.0e-3 * currentDualBound_);
1053                if (upperSaved_[iSequence] - abcUpper[iSequence] > abcLower[iSequence] - lowerSaved_[iSequence]) {
1054                  setFakeBound(iSequence, AbcSimplexDual::upperFake);
1055                  //assert (abcLower[iSequence]==lowerSaved_[iSequence]);
1056                  abcLower[iSequence] = CoinMax(abcLower[iSequence], lowerSaved_[iSequence]);
1057                  assert(abcUpper[iSequence] < upperSaved_[iSequence]);
1058                } else {
1059                  setFakeBound(iSequence, AbcSimplexDual::lowerFake);
1060                  assert(abcLower[iSequence] > lowerSaved_[iSequence]);
1061                  //assert (abcUpper[iSequence]==upperSaved_[iSequence]);
1062                  abcUpper[iSequence] = CoinMin(abcUpper[iSequence], upperSaved_[iSequence]);
1063                }
1064              }
1065            }
1066          } else {
1067            if (getFakeBound(iSequence) == AbcSimplexDual::noFake)
1068              numberFake_++;
1069            abcLower[iSequence] = -0.5 * currentDualBound_;
1070            abcUpper[iSequence] = 0.5 * currentDualBound_;
1071            setFakeBound(iSequence, AbcSimplexDual::bothFake);
1072            abort();
1073          }
1074          if (status == atUpperBound)
1075            abcSolution[iSequence] = abcUpper[iSequence];
1076          else
1077            abcSolution[iSequence] = abcLower[iSequence];
1078        } else {
1079          // set non basic free variables to fake bounds
1080          // I don't think we should ever get here
1081          CoinAssert(!("should not be here"));
1082          abcLower[iSequence] = -0.5 * currentDualBound_;
1083          abcUpper[iSequence] = 0.5 * currentDualBound_;
1084          setFakeBound(iSequence, AbcSimplexDual::bothFake);
1085          numberFake_++;
1086          setInternalStatus(iSequence, atUpperBound);
1087          abcSolution[iSequence] = 0.5 * currentDualBound_;
1088        }
1089      } else if (status == basic) {
1090        // make sure not at fake bound and bounds correct
1091        setFakeBound(iSequence, AbcSimplexDual::noFake);
1092        double gap = abcUpper[iSequence] - abcLower[iSequence];
1093        if (gap > 0.5 * currentDualBound_ && gap < 2.0 * currentDualBound_) {
1094          abcLower[iSequence] = lowerSaved_[iSequence];
1095          ;
1096          abcUpper[iSequence] = upperSaved_[iSequence];
1097          ;
1098        }
1099      }
1100    }
1101    return 0;
1102  }
1103}
1104/*
1105  order of variables is a) atLo/atUp
1106                        b) free/superBasic
1107                        pre pass done in AbcMatrix to get possibles
1108                        dj's will be stored in a parallel array (tempArray) for first pass
1109                        this pass to be done inside AbcMatrix then as now
1110                        use perturbationArray as variable tolerance
1111  idea is to use large tolerances and then balance and reduce every so often
1112  this should be in conjunction with dualBound changes
1113  variables always stay at bounds on accuracy issues - only deliberate flips allowed
1114  variable tolerance should be 0.8-1.0 desired (at end)
1115  use going to clpMatrix_ as last resort
1116 */
1117/*
1118  Array has tableau row
1119  Puts candidates in list
1120  Returns upper theta (infinity if no pivot) and may set sequenceIn_ if free
1121*/
1122void AbcSimplexDual::dualColumn1(bool doAll)
1123{
1124  const CoinIndexedVector &update = usefulArray_[arrayForBtran_];
1125  CoinPartitionedVector &tableauRow = usefulArray_[arrayForTableauRow_];
1126  CoinPartitionedVector &candidateList = usefulArray_[arrayForDualColumn_];
1127  if ((stateOfProblem_ & VALUES_PASS) == 0)
1128    upperTheta_ = 1.0e31;
1129  else
1130    upperTheta_ = fabs(abcDj_[sequenceOut_]) - 0.5 * dualTolerance_;
1131  assert(upperTheta_ > 0.0);
1132  if (!doAll)
1133    upperTheta_ = abcMatrix_->dualColumn1(update, tableauRow, candidateList);
1134  else
1135    upperTheta_ = dualColumn1B();
1136  //tableauRow.compact();
1137  //candidateList.compact();
1138  //tableauRow.setPackedMode(true);
1139  //candidateList.setPackedMode(true);
1140}
1141/*
1142  Array has tableau row
1143  Puts candidates in list
1144  Returns upper theta (infinity if no pivot) and may set sequenceIn_ if free
1145*/
1146double
1147AbcSimplexDual::dualColumn1A()
1148{
1149  const CoinIndexedVector &update = usefulArray_[arrayForBtran_];
1150  CoinPartitionedVector &tableauRow = usefulArray_[arrayForTableauRow_];
1151  CoinPartitionedVector &candidateList = usefulArray_[arrayForDualColumn_];
1152  int numberRemaining = 0;
1153  double upperTheta = upperTheta_;
1154  // pivot elements
1155  double *COIN_RESTRICT arrayCandidate = candidateList.denseVector();
1156  // indices
1157  int *COIN_RESTRICT indexCandidate = candidateList.getIndices();
1158  const double *COIN_RESTRICT abcDj = abcDj_;
1159  //const double *  COIN_RESTRICT abcPerturbation = abcPerturbation_;
1160  const unsigned char *COIN_RESTRICT internalStatus = internalStatus_;
1161  const double *COIN_RESTRICT pi = update.denseVector();
1162  const int *COIN_RESTRICT piIndex = update.getIndices();
1163  int number = update.getNumElements();
1164  int *COIN_RESTRICT index = tableauRow.getIndices();
1165  double *COIN_RESTRICT array = tableauRow.denseVector();
1166  assert(!tableauRow.getNumElements());
1167  int numberNonZero = 0;
1168  // do first pass to get possibles
1169  double bestPossible = 0.0;
1170  // We can also see if infeasible or pivoting on free
1171  double tentativeTheta = 1.0e25; // try with this much smaller as guess
1172  double acceptablePivot = currentAcceptablePivot_;
1173  double dualT = -currentDualTolerance_;
1174  const double multiplier[] = { 1.0, -1.0 };
1175  if (ordinaryVariables_) {
1176    for (int i = 0; i < number; i++) {
1177      int iRow = piIndex[i];
1178      unsigned char type = internalStatus[iRow];
1179      if ((type & 4) == 0) {
1180        index[numberNonZero] = iRow;
1181        double piValue = pi[iRow];
1182        array[numberNonZero++] = piValue;
1183        int iStatus = internalStatus[iRow] & 3;
1184        double mult = multiplier[iStatus];
1185        double alpha = piValue * mult;
1186        double oldValue = abcDj[iRow] * mult;
1187        double value = oldValue - tentativeTheta * alpha;
1188        if (value < dualT) {
1189          bestPossible = CoinMax(bestPossible, alpha);
1190          value = oldValue - upperTheta * alpha;
1191          if (value < dualT && alpha >= acceptablePivot) {
1192            upperTheta = (oldValue - dualT) / alpha;
1193          }
1194          // add to list
1195          arrayCandidate[numberRemaining] = alpha;
1196          indexCandidate[numberRemaining++] = iRow;
1197        }
1198      }
1199    }
1200  } else {
1201    double badFree = 0.0;
1202    double freePivot = currentAcceptablePivot_;
1203    alpha_ = 0.0;
1204    // We can also see if infeasible or pivoting on free
1205    for (int i = 0; i < number; i++) {
1206      int iRow = piIndex[i];
1207      unsigned char type = internalStatus[iRow];
1208      double piValue = pi[iRow];
1209      if ((type & 4) == 0) {
1210        index[numberNonZero] = iRow;
1211        array[numberNonZero++] = piValue;
1212        int iStatus = internalStatus[iRow] & 3;
1213        double mult = multiplier[iStatus];
1214        double alpha = piValue * mult;
1215        double oldValue = abcDj[iRow] * mult;
1216        double value = oldValue - tentativeTheta * alpha;
1217        if (value < dualT) {
1218          bestPossible = CoinMax(bestPossible, alpha);
1219          value = oldValue - upperTheta * alpha;
1220          if (value < dualT && alpha >= acceptablePivot) {
1221            upperTheta = (oldValue - dualT) / alpha;
1222          }
1223          // add to list
1224          arrayCandidate[numberRemaining] = alpha;
1225          indexCandidate[numberRemaining++] = iRow;
1226        }
1227      } else if ((type & 7) < 6) {
1228        //} else if (getInternalStatus(iRow)==isFree||getInternalStatus(iRow)==superBasic) {
1229        bool keep;
1230        index[numberNonZero] = iRow;
1231        array[numberNonZero++] = piValue;
1232        bestPossible = CoinMax(bestPossible, fabs(piValue));
1233        double oldValue = abcDj[iRow];
1234        // If free has to be very large - should come in via dualRow
1235        //if (getInternalStatus(iRow+addSequence)==isFree&&fabs(piValue)<1.0e-3)
1236        //break;
1237        if (oldValue > currentDualTolerance_) {
1238          keep = true;
1239        } else if (oldValue < -currentDualTolerance_) {
1240          keep = true;
1241        } else {
1242          if (fabs(piValue) > CoinMax(10.0 * currentAcceptablePivot_, 1.0e-5)) {
1243            keep = true;
1244          } else {
1245            keep = false;
1246            badFree = CoinMax(badFree, fabs(piValue));
1247          }
1248        }
1249        if (keep) {
1250          // free - choose largest
1251          if (fabs(piValue) > freePivot) {
1252            freePivot = fabs(piValue);
1253            sequenceIn_ = iRow;
1254            theta_ = oldValue / piValue;
1255            alpha_ = piValue;
1256          }
1257        }
1258      }
1259    }
1260  }
1261  //tableauRow.setNumElements(numberNonZero);
1262  //candidateList.setNumElements(numberRemaining);
1263  tableauRow.setNumElementsPartition(0, numberNonZero);
1264  candidateList.setNumElementsPartition(0, numberRemaining);
1265  if (!numberRemaining && sequenceIn_ < 0 && upperTheta_ > 1.0e29) {
1266    return COIN_DBL_MAX; // Looks infeasible
1267  } else {
1268    return upperTheta;
1269  }
1270}
1271/*
1272  Array has tableau row
1273  Puts candidates in list
1274  Returns upper theta (infinity if no pivot) and may set sequenceIn_ if free
1275*/
1276double
1277AbcSimplexDual::dualColumn1B()
1278{
1279  CoinPartitionedVector &tableauRow = usefulArray_[arrayForTableauRow_];
1280  CoinPartitionedVector &candidateList = usefulArray_[arrayForDualColumn_];
1281  tableauRow.compact();
1282  candidateList.clearAndReset();
1283  tableauRow.setPackedMode(true);
1284  candidateList.setPackedMode(true);
1285  int numberRemaining = 0;
1286  double upperTheta = upperTheta_;
1287  // pivot elements
1288  double *COIN_RESTRICT arrayCandidate = candidateList.denseVector();
1289  // indices
1290  int *COIN_RESTRICT indexCandidate = candidateList.getIndices();
1291  const double *COIN_RESTRICT abcDj = abcDj_;
1292  const unsigned char *COIN_RESTRICT internalStatus = internalStatus_;
1293  int *COIN_RESTRICT index = tableauRow.getIndices();
1294  double *COIN_RESTRICT array = tableauRow.denseVector();
1295  int numberNonZero = tableauRow.getNumElements();
1296  // do first pass to get possibles
1297  double bestPossible = 0.0;
1298  // We can also see if infeasible or pivoting on free
1299  double tentativeTheta = 1.0e25; // try with this much smaller as guess
1300  double acceptablePivot = currentAcceptablePivot_;
1301  double dualT = -currentDualTolerance_;
1302  const double multiplier[] = { 1.0, -1.0 };
1303  if (ordinaryVariables_) {
1304    for (int i = 0; i < numberNonZero; i++) {
1305      int iSequence = index[i];
1306      unsigned char iStatus = internalStatus[iSequence] & 7;
1307      if (iStatus < 4) {
1308        double tableauValue = array[i];
1309        if (fabs(tableauValue) > zeroTolerance_) {
1310          double mult = multiplier[iStatus];
1311          double alpha = tableauValue * mult;
1312          double oldValue = abcDj[iSequence] * mult;
1313          double value = oldValue - tentativeTheta * alpha;
1314          if (value < dualT) {
1315            bestPossible = CoinMax(bestPossible, alpha);
1316            value = oldValue - upperTheta * alpha;
1317            if (value < dualT && alpha >= acceptablePivot) {
1318              upperTheta = (oldValue - dualT) / alpha;
1319            }
1320            // add to list
1321            arrayCandidate[numberRemaining] = alpha;
1322            indexCandidate[numberRemaining++] = iSequence;
1323          }
1324        }
1325      }
1326    }
1327  } else {
1328    double badFree = 0.0;
1329    double freePivot = currentAcceptablePivot_;
1330    double currentDualTolerance = currentDualTolerance_;
1331    for (int i = 0; i < numberNonZero; i++) {
1332      int iSequence = index[i];
1333      unsigned char iStatus = internalStatus[iSequence] & 7;
1334      if (iStatus < 6) {
1335        double tableauValue = array[i];
1336        if (fabs(tableauValue) > zeroTolerance_) {
1337          if ((iStatus & 4) == 0) {
1338            double mult = multiplier[iStatus];
1339            double alpha = tableauValue * mult;
1340            double oldValue = abcDj[iSequence] * mult;
1341            double value = oldValue - tentativeTheta * alpha;
1342            if (value < dualT) {
1343              bestPossible = CoinMax(bestPossible, alpha);
1344              value = oldValue - upperTheta * alpha;
1345              if (value < dualT && alpha >= acceptablePivot) {
1346                upperTheta = (oldValue - dualT) / alpha;
1347              }
1348              // add to list
1349              arrayCandidate[numberRemaining] = alpha;
1350              indexCandidate[numberRemaining++] = iSequence;
1351            }
1352          } else {
1353            bool keep;
1354            bestPossible = CoinMax(bestPossible, fabs(tableauValue));
1355            double oldValue = abcDj[iSequence];
1356            // If free has to be very large - should come in via dualRow
1357            //if (getInternalStatus(iSequence+addSequence)==isFree&&fabs(tableauValue)<1.0e-3)
1358            //break;
1359            if (oldValue > currentDualTolerance) {
1360              keep = true;
1361            } else if (oldValue < -currentDualTolerance) {
1362              keep = true;
1363            } else {
1364              if (fabs(tableauValue) > CoinMax(10.0 * acceptablePivot, 1.0e-5)) {
1365                keep = true;
1366              } else {
1367                keep = false;
1368                badFree = CoinMax(badFree, fabs(tableauValue));
1369              }
1370            }
1371            if (keep) {
1372#ifdef PAN
1373              if (fakeSuperBasic(iSequence) >= 0) {
1374#endif
1375                if (iSequence == freeSequenceIn_)
1376                  tableauValue = COIN_DBL_MAX;
1377                // free - choose largest
1378                if (fabs(tableauValue) > fabs(freePivot)) {
1379                  freePivot = tableauValue;
1380                  sequenceIn_ = iSequence;
1381                  theta_ = oldValue / tableauValue;
1382                  alpha_ = tableauValue;
1383                }
1384#ifdef PAN
1385              }
1386#endif
1387            }
1388          }
1389        }
1390      }
1391    }
1392  }
1393  candidateList.setNumElements(numberRemaining);
1394  if (!numberRemaining && sequenceIn_ < 0 && upperTheta_ > 1.0e29) {
1395    return COIN_DBL_MAX; // Looks infeasible
1396  } else {
1397    return upperTheta;
1398  }
1399}
1400//#define MORE_CAREFUL
1401//#define PRINT_RATIO_PROGRESS
1402#define MINIMUMTHETA 1.0e-18
1403#define MAXTRY 100
1404#define TOL_TYPE -1
1405#if TOL_TYPE == -1
1406#define useTolerance1() (0.0)
1407#define useTolerance2() (0.0)
1408#define goToTolerance1() (0.0)
1409#define goToTolerance2() (0.0)
1410#elif TOL_TYPE == 0
1411#define useTolerance1() (currentDualTolerance_)
1412#define useTolerance2() (currentDualTolerance_)
1413#define goToTolerance1() (0.0)
1414#define goToTolerance2() (0.0)
1415#elif TOL_TYPE == 1
1416#define useTolerance1() (perturbedTolerance)
1417#define useTolerance2() (0.0)
1418#define goToTolerance1() (0.0)
1419#define goToTolerance2() (0.0)
1420#elif TOL_TYPE == 2
1421#define useTolerance1() (perturbedTolerance)
1422#define useTolerance2() (perturbedTolerance)
1423#define goToTolerance1() (0.0)
1424#define goToTolerance2() (0.0)
1425#elif TOL_TYPE == 3
1426#define useTolerance1() (perturbedTolerance)
1427#define useTolerance2() (perturbedTolerance)
1428#define goToTolerance1() (perturbedTolerance)
1429#define goToTolerance2() (0.0)
1430#elif TOL_TYPE == 4
1431#define useTolerance1() (perturbedTolerance)
1432#define useTolerance2() (perturbedTolerance)
1433#define goToTolerance1() (perturbedTolerance)
1434#define goToTolerance2() (perturbedTolerance)
1435#else
1436XXXX
1437#endif
1438#if CLP_CAN_HAVE_ZERO_OBJ < 2
1439#define MODIFYCOST 2
1440#endif
1441#define DONT_MOVE_OBJECTIVE
1442static void dualColumn2Bit(AbcSimplexDual *dual, dualColumnResult *result,
1443  int numberBlocks)
1444{
1445  for (int iBlock = 1; iBlock < numberBlocks; iBlock++) {
1446    cilk_spawn dual->dualColumn2First(result[iBlock]);
1447  }
1448  dual->dualColumn2First(result[0]);
1449  cilk_sync;
1450}
1451void AbcSimplexDual::dualColumn2First(dualColumnResult &result)
1452{
1453  CoinPartitionedVector &candidateList = usefulArray_[arrayForDualColumn_];
1454  CoinPartitionedVector &flipList = usefulArray_[arrayForFlipBounds_];
1455  int iBlock = result.block;
1456  int numberSwapped = result.numberSwapped;
1457  int numberRemaining = result.numberRemaining;
1458  double tentativeTheta = result.tentativeTheta;
1459  // pivot elements
1460  double *COIN_RESTRICT array = candidateList.denseVector();
1461  // indices
1462  int *COIN_RESTRICT indices = candidateList.getIndices();
1463  const double *COIN_RESTRICT abcLower = abcLower_;
1464  const double *COIN_RESTRICT abcUpper = abcUpper_;
1465  const double *COIN_RESTRICT abcDj = abcDj_;
1466  const unsigned char *COIN_RESTRICT internalStatus = internalStatus_;
1467  const double multiplier[] = { 1.0, -1.0 };
1468
1469  // For list of flipped
1470  int *COIN_RESTRICT flippedIndex = flipList.getIndices();
1471  double *COIN_RESTRICT flippedElement = flipList.denseVector();
1472  // adjust
1473  //if (iBlock>=0) {
1474  int offset = candidateList.startPartition(iBlock);
1475  assert(offset == flipList.startPartition(iBlock));
1476  array += offset;
1477  indices += offset;
1478  flippedIndex += offset;
1479  flippedElement += offset;
1480  //}
1481  double thruThis = 0.0;
1482
1483  double bestPivot = currentAcceptablePivot_;
1484  int bestSequence = -1;
1485  int numberInList = numberRemaining;
1486  numberRemaining = 0;
1487  double upperTheta = 1.0e50;
1488  double increaseInThis = 0.0; //objective increase in this loop
1489  for (int i = 0; i < numberInList; i++) {
1490    int iSequence = indices[i];
1491    assert(getInternalStatus(iSequence) != isFree
1492      && getInternalStatus(iSequence) != superBasic);
1493    int iStatus = internalStatus[iSequence] & 3;
1494    // treat as if at lower bound
1495    double mult = multiplier[iStatus];
1496    double alpha = array[i];
1497    double oldValue = abcDj[iSequence] * mult;
1498    double value = oldValue - tentativeTheta * alpha;
1499    assert(alpha > 0.0);
1500    //double perturbedTolerance = abcPerturbation[iSequence];
1501    if (value < -currentDualTolerance_) {
1502      // possible here or could go through
1503      double range = abcUpper[iSequence] - abcLower[iSequence];
1504      thruThis += range * alpha;
1505      increaseInThis += (oldValue - useTolerance1()) * range;
1506      // goes on swapped list (also means candidates if too many)
1507      flippedElement[numberSwapped] = alpha;
1508      flippedIndex[numberSwapped++] = iSequence;
1509      if (alpha > bestPivot) {
1510        bestPivot = alpha;
1511        bestSequence = iSequence;
1512      }
1513    } else {
1514      // won't go through this time
1515      value = oldValue - upperTheta * alpha;
1516      if (value < -currentDualTolerance_ && alpha >= currentAcceptablePivot_)
1517        upperTheta = (oldValue + currentDualTolerance_) / alpha;
1518      array[numberRemaining] = alpha;
1519      indices[numberRemaining++] = iSequence;
1520    }
1521  }
1522  result.bestEverPivot = bestPivot;
1523  result.sequence = bestSequence;
1524  result.numberSwapped = numberSwapped;
1525  result.numberRemaining = numberRemaining;
1526  result.thruThis = thruThis;
1527  result.increaseInThis = increaseInThis;
1528  result.theta = upperTheta;
1529}
1530void AbcSimplexDual::dualColumn2Most(dualColumnResult &result)
1531{
1532  CoinPartitionedVector &candidateList = usefulArray_[arrayForDualColumn_];
1533  CoinPartitionedVector &flipList = usefulArray_[arrayForFlipBounds_];
1534  int numberBlocksReal = candidateList.getNumPartitions();
1535  flipList.setPartitions(numberBlocksReal, candidateList.startPartitions());
1536  int numberBlocks;
1537  if (!numberBlocksReal)
1538    numberBlocks = 1;
1539  else
1540    numberBlocks = numberBlocksReal;
1541  //usefulArray_[arrayForTableauRow_].compact();
1542  //candidateList.compact();
1543  int numberSwapped = 0;
1544  int numberRemaining = candidateList.getNumElements();
1545
1546  double totalThru = 0.0; // for when variables flip
1547  int sequenceIn = -1;
1548  double bestEverPivot = currentAcceptablePivot_;
1549  // If we think we need to modify costs (not if something from broad sweep)
1550  bool modifyCosts = false;
1551  // Increase in objective due to swapping bounds (may be negative)
1552  double increaseInObjective = 0.0;
1553  // pivot elements
1554  double *COIN_RESTRICT array = candidateList.denseVector();
1555  // indices
1556  int *COIN_RESTRICT indices = candidateList.getIndices();
1557  const double *COIN_RESTRICT abcLower = abcLower_;
1558  const double *COIN_RESTRICT abcUpper = abcUpper_;
1559  //const double *  COIN_RESTRICT abcSolution = abcSolution_;
1560  const double *COIN_RESTRICT abcDj = abcDj_;
1561  //const double *  COIN_RESTRICT abcPerturbation = abcPerturbation_;
1562  const unsigned char *COIN_RESTRICT internalStatus = internalStatus_;
1563  // We can also see if infeasible or pivoting on free
1564  double tentativeTheta = 1.0e25; // try with this much smaller as guess
1565  const double multiplier[] = { 1.0, -1.0 };
1566
1567  // For list of flipped
1568  int numberLastSwapped = 0;
1569  int *COIN_RESTRICT flippedIndex = flipList.getIndices();
1570  double *COIN_RESTRICT flippedElement = flipList.denseVector();
1571  // If sum of bad small pivots too much
1572#ifdef MORE_CAREFUL
1573  bool badSumPivots = false;
1574#endif
1575  // not free
1576  int lastSequence = -1;
1577  double lastTheta = 0.0;
1578  double lastPivotValue = 0.0;
1579  double thisPivotValue = 0.0;
1580  // if free then always choose , otherwise ...
1581  double theta = 1.0e50;
1582  // now flip flop between spare arrays until reasonable theta
1583  tentativeTheta = CoinMax(10.0 * upperTheta_, 1.0e-7);
1584
1585  // loops increasing tentative theta until can't go through
1586#ifdef PRINT_RATIO_PROGRESS
1587  int iPass = 0;
1588#endif
1589  dualColumnResult result2[8];
1590  for (int iBlock = 0; iBlock < numberBlocks; iBlock++) {
1591    result2[iBlock].block = iBlock;
1592    if (numberBlocksReal)
1593      result2[iBlock].numberRemaining = candidateList.getNumElements(iBlock);
1594    else
1595      result2[iBlock].numberRemaining = candidateList.getNumElements();
1596    result2[iBlock].numberSwapped = 0;
1597    result2[iBlock].numberLastSwapped = 0;
1598  }
1599  double useThru = 0.0;
1600  while (tentativeTheta < 1.0e22) {
1601    double thruThis = 0.0;
1602
1603    double bestPivot = currentAcceptablePivot_;
1604    int bestSequence = -1;
1605    sequenceIn = -1;
1606    thisPivotValue = 0.0;
1607    int numberInList = numberRemaining;
1608    double upperTheta = 1.0e50;
1609    double increaseInThis = 0.0; //objective increase in this loop
1610    for (int iBlock = 0; iBlock < numberBlocks; iBlock++) {
1611      result2[iBlock].tentativeTheta = tentativeTheta;
1612      //result2[iBlock].numberRemaining=numberRemaining;
1613      //result2[iBlock].numberSwapped=numberSwapped;
1614      //result2[iBlock].numberLastSwapped=numberLastSwapped;
1615    }
1616    for (int iBlock = 1; iBlock < numberBlocks; iBlock++) {
1617      cilk_spawn dualColumn2First(result2[iBlock]);
1618    }
1619    dualColumn2First(result2[0]);
1620    cilk_sync;
1621    //dualColumn2Bit(this,result2,numberBlocks);
1622    numberSwapped = 0;
1623    numberRemaining = 0;
1624    for (int iBlock = 0; iBlock < numberBlocks; iBlock++) {
1625      if (result2[iBlock].bestEverPivot > bestPivot) {
1626        bestPivot = result2[iBlock].bestEverPivot;
1627        bestSequence = result2[iBlock].sequence;
1628      }
1629      numberSwapped += result2[iBlock].numberSwapped;
1630      numberRemaining += result2[iBlock].numberRemaining;
1631      thruThis += result2[iBlock].thruThis;
1632      increaseInThis += result2[iBlock].increaseInThis;
1633      if (upperTheta > result2[iBlock].theta)
1634        upperTheta = result2[iBlock].theta;
1635    }
1636    // end of pass
1637#ifdef PRINT_RATIO_PROGRESS
1638    iPass++;
1639    printf("Iteration %d pass %d dualOut %g thru %g increase %g numberRemain %d\n",
1640      numberIterations_, iPass, fabs(dualOut_), totalThru + thruThis, increaseInObjective + increaseInThis,
1641      numberRemaining);
1642#endif
1643    if (totalThru + thruThis > fabs(dualOut_) - fabs(1.0e-9 * dualOut_) - 1.0e-12 || (increaseInObjective + increaseInThis < 0.0 && bestSequence >= 0) || !numberRemaining) {
1644      // We should be pivoting in this batch
1645      // so compress down to this lot
1646      numberRemaining = 0;
1647      candidateList.clearAndReset();
1648      for (int iBlock = 0; iBlock < numberBlocks; iBlock++) {
1649        int offset = flipList.startPartition(iBlock);
1650        // set so can compact
1651        int numberKeep = result2[iBlock].numberLastSwapped;
1652        flipList.setNumElementsPartition(iBlock, numberKeep);
1653        for (int i = offset + numberKeep; i < offset + result2[iBlock].numberSwapped; i++) {
1654          array[numberRemaining] = flippedElement[i];
1655          flippedElement[i] = 0;
1656          indices[numberRemaining++] = flippedIndex[i];
1657        }
1658      }
1659      // make sure candidateList will be cleared correctly
1660      candidateList.setNumElements(numberRemaining);
1661      candidateList.setPackedMode(true);
1662      flipList.compact();
1663      numberSwapped = numberLastSwapped;
1664      useThru = totalThru;
1665      assert(totalThru <= fabs(dualOut_));
1666      int iTry;
1667      // first get ratio with tolerance
1668      for (iTry = 0; iTry < MAXTRY; iTry++) {
1669
1670        upperTheta = 1.0e50;
1671        numberInList = numberRemaining;
1672        numberRemaining = 0;
1673        increaseInThis = 0.0; //objective increase in this loop
1674        thruThis = 0.0;
1675        for (int i = 0; i < numberInList; i++) {
1676          int iSequence = indices[i];
1677          assert(getInternalStatus(iSequence) != isFree
1678            && getInternalStatus(iSequence) != superBasic);
1679          int iStatus = internalStatus[iSequence] & 3;
1680          // treat as if at lower bound
1681          double mult = multiplier[iStatus];
1682          double alpha = array[i];
1683          double oldValue = abcDj[iSequence] * mult;
1684          double value = oldValue - upperTheta * alpha;
1685          assert(alpha > 0.0);
1686          //double perturbedTolerance = abcPerturbation[iSequence];
1687          if (value < -currentDualTolerance_) {
1688            if (alpha >= currentAcceptablePivot_) {
1689              upperTheta = (oldValue + currentDualTolerance_) / alpha;
1690            }
1691          }
1692        }
1693        bestPivot = currentAcceptablePivot_;
1694        sequenceIn = -1;
1695        double largestPivot = currentAcceptablePivot_;
1696        // Sum of bad small pivots
1697#ifdef MORE_CAREFUL
1698        double sumBadPivots = 0.0;
1699        badSumPivots = false;
1700#endif
1701        // Make sure upperTheta will work (-O2 and above gives problems)
1702        upperTheta *= 1.0000000001;
1703        for (int i = 0; i < numberInList; i++) {
1704          int iSequence = indices[i];
1705          int iStatus = internalStatus[iSequence] & 3;
1706          // treat as if at lower bound
1707          double mult = multiplier[iStatus];
1708          double alpha = array[i];
1709          double oldValue = abcDj[iSequence] * mult;
1710          double value = oldValue - upperTheta * alpha;
1711          assert(alpha > 0.0);
1712          //double perturbedTolerance = abcPerturbation[iSequence];
1713          double badDj = 0.0;
1714          if (value <= 0.0) {
1715            // add to list of swapped
1716            flippedElement[numberSwapped] = alpha;
1717            flippedIndex[numberSwapped++] = iSequence;
1718            badDj = oldValue - useTolerance2();
1719            // select if largest pivot
1720            bool take = (alpha > bestPivot);
1721#ifdef MORE_CAREFUL
1722            if (alpha < currentAcceptablePivot_ && upperTheta < 1.0e20) {
1723              if (value < -currentDualTolerance_) {
1724                double gap = abcUpper[iSequence] - abcLower[iSequence];
1725                if (gap < 1.0e20)
1726                  sumBadPivots -= value * gap;
1727                else
1728                  sumBadPivots += 1.0e20;
1729                //printf("bad %d alpha %g dj at lower %g\n",
1730                //     iSequence,alpha,value);
1731              }
1732            }
1733#endif
1734            if (take && flagged(iSequence)) {
1735              if (bestPivot > currentAcceptablePivot_)
1736                take = false;
1737            }
1738            if (take) {
1739              sequenceIn = iSequence;
1740              thisPivotValue = alpha;
1741              bestPivot = alpha;
1742              if (flagged(iSequence))
1743                bestPivot = currentAcceptablePivot_;
1744              bestSequence = iSequence;
1745              theta = (oldValue + goToTolerance1()) / alpha;
1746              //assert (oldValue==abcDj[iSequence]);
1747              largestPivot = CoinMax(largestPivot, 0.5 * bestPivot);
1748            }
1749            double range = abcUpper[iSequence] - abcLower[iSequence];
1750            thruThis += range * fabs(alpha);
1751            increaseInThis += badDj * range;
1752          } else {
1753            // add to list of remaining
1754            array[numberRemaining] = alpha;
1755            indices[numberRemaining++] = iSequence;
1756          }
1757        }
1758#ifdef MORE_CAREFUL
1759        // If we have done pivots and things look bad set alpha_ 0.0 to force factorization
1760        if (sumBadPivots > 1.0e4) {
1761#if ABC_NORMAL_DEBUG > 0
1762          if (handler_->logLevel() > 1)
1763            printf("maybe forcing re-factorization - sum %g  %d pivots\n", sumBadPivots,
1764              abcFactorization_->pivots());
1765#endif
1766          if (abcFactorization_->pivots() > 3) {
1767            badSumPivots = true;
1768#ifdef PRINT_RATIO_PROGRESS
1769            printf("maybe forcing re-factorization - sum %g  %d pivots\n", sumBadPivots,
1770              abcFactorization_->pivots());
1771#endif
1772            break;
1773          }
1774        }
1775#endif
1776        // If we stop now this will be increase in objective (I think)
1777        double increase = (fabs(dualOut_) - totalThru) * theta;
1778        increase += increaseInObjective;
1779        if (theta < 0.0)
1780          thruThis += fabs(dualOut_); // force using this one
1781        if (((increaseInObjective < 0.0 && increase < 0.0) || (theta > 1.0e9 && lastTheta < 1.0e-5))
1782          && lastSequence >= 0) {
1783          // back
1784          // We may need to be more careful - we could do by
1785          // switch so we always do fine grained?
1786          bestPivot = 0.0;
1787          bestSequence = -1;
1788        } else {
1789          // add in
1790          totalThru += thruThis;
1791          increaseInObjective += increaseInThis;
1792          lastTheta = theta;
1793        }
1794        if (bestPivot < 0.1 * bestEverPivot && bestEverPivot > 1.0e-6 && (bestPivot < 1.0e-3 || totalThru * 2.0 > fabs(dualOut_))) {
1795          // back to previous one
1796          sequenceIn = lastSequence;
1797          thisPivotValue = lastPivotValue;
1798#ifdef PRINT_RATIO_PROGRESS
1799          printf("Try %d back to previous bestPivot %g bestEver %g totalThru %g\n",
1800            iTry, bestPivot, bestEverPivot, totalThru);
1801#endif
1802          break;
1803        } else if (sequenceIn == -1 && upperTheta > largeValue_) {
1804          if (lastPivotValue > currentAcceptablePivot_) {
1805            // back to previous one
1806            sequenceIn = lastSequence;
1807            thisPivotValue = lastPivotValue;
1808#ifdef PRINT_RATIO_PROGRESS
1809            printf("Try %d no pivot this time large theta %g lastPivot %g\n",
1810              iTry, upperTheta, lastPivotValue);
1811#endif
1812          } else {
1813            // can only get here if all pivots too small
1814#ifdef PRINT_RATIO_PROGRESS
1815            printf("XXBAD Try %d no pivot this time large theta %g lastPivot %g\n",
1816              iTry, upperTheta, lastPivotValue);
1817#endif
1818          }
1819          break;
1820        } else if (totalThru >= fabs(dualOut_)) {
1821#ifdef PRINT_RATIO_PROGRESS
1822          printf("Try %d upperTheta %g totalThru %g - can modify costs\n",
1823            iTry, upperTheta, lastPivotValue);
1824#endif
1825          modifyCosts = true; // fine grain - we can modify costs
1826          break; // no point trying another loop
1827        } else {
1828          lastSequence = sequenceIn;
1829          lastPivotValue = thisPivotValue;
1830          if (bestPivot > bestEverPivot) {
1831            bestEverPivot = bestPivot;
1832            //bestEverSequence=bestSequence;
1833          }
1834          modifyCosts = true; // fine grain - we can modify costs
1835        }
1836#ifdef PRINT_RATIO_PROGRESS
1837        printf("Onto next pass totalthru %g bestPivot %g\n",
1838          totalThru, lastPivotValue);
1839#endif
1840        numberLastSwapped = numberSwapped;
1841      }
1842      break;
1843    } else {
1844      // skip this lot
1845      if (bestPivot > 1.0e-3 || bestPivot > bestEverPivot) {
1846#ifdef PRINT_RATIO_PROGRESS
1847        printf("Continuing bestPivot %g bestEver %g %d swapped - new tentative %g\n",
1848          bestPivot, bestEverPivot, numberSwapped, 2 * upperTheta);
1849#endif
1850        bestEverPivot = bestPivot;
1851        //bestEverSequence=bestSequence;
1852        lastPivotValue = bestPivot;
1853        lastSequence = bestSequence;
1854        for (int iBlock = 0; iBlock < numberBlocks; iBlock++)
1855          result2[iBlock].numberLastSwapped = result2[iBlock].numberSwapped;
1856        numberLastSwapped = numberSwapped;
1857      } else {
1858#ifdef PRINT_RATIO_PROGRESS
1859        printf("keep old swapped bestPivot %g bestEver %g %d swapped - new tentative %g\n",
1860          bestPivot, bestEverPivot, numberSwapped, 2 * upperTheta);
1861#endif
1862        // keep old swapped
1863        numberRemaining = 0;
1864        for (int iBlock = 0; iBlock < numberBlocks; iBlock++) {
1865          int offset = flipList.startPartition(iBlock);
1866          int numberRemaining2 = offset + result2[iBlock].numberRemaining;
1867          int numberLastSwapped = result2[iBlock].numberLastSwapped;
1868          for (int i = offset + numberLastSwapped;
1869               i < offset + result2[iBlock].numberSwapped; i++) {
1870            array[numberRemaining2] = flippedElement[i];
1871            flippedElement[i] = 0;
1872            indices[numberRemaining2++] = flippedIndex[i];
1873          }
1874          result2[iBlock].numberRemaining = numberRemaining2 - offset;
1875          numberRemaining += result2[iBlock].numberRemaining;
1876          result2[iBlock].numberSwapped = numberLastSwapped;
1877        }
1878        numberSwapped = numberLastSwapped;
1879      }
1880      increaseInObjective += increaseInThis;
1881      tentativeTheta = 2.0 * upperTheta;
1882      totalThru += thruThis;
1883    }
1884  }
1885  if (flipList.getNumPartitions()) {
1886    // need to compress
1887    numberRemaining = 0;
1888    candidateList.clearAndReset();
1889    for (int iBlock = 0; iBlock < numberBlocks; iBlock++) {
1890      int offset = flipList.startPartition(iBlock);
1891      // set so can compact
1892      int numberKeep = result2[iBlock].numberLastSwapped;
1893      flipList.setNumElementsPartition(iBlock, numberKeep);
1894      for (int i = offset + numberKeep; i < offset + result2[iBlock].numberSwapped; i++) {
1895        array[numberRemaining] = flippedElement[i];
1896        flippedElement[i] = 0;
1897        indices[numberRemaining++] = flippedIndex[i];
1898      }
1899    }
1900    // make sure candidateList will be cleared correctly
1901    candidateList.setNumElements(numberRemaining);
1902    candidateList.setPackedMode(true);
1903    flipList.compact();
1904    numberSwapped = numberLastSwapped;
1905  }
1906  result.theta = theta; //?
1907  result.totalThru = totalThru;
1908  result.useThru = useThru;
1909  result.bestEverPivot = bestEverPivot; //?
1910  result.increaseInObjective = increaseInObjective; //?
1911  result.tentativeTheta = tentativeTheta; //?
1912  result.lastPivotValue = lastPivotValue;
1913  result.thisPivotValue = thisPivotValue;
1914  result.lastSequence = lastSequence;
1915  result.sequence = sequenceIn;
1916  //result.block;
1917  //result.pass;
1918  //result.phase;
1919  result.numberSwapped = numberSwapped;
1920  result.numberLastSwapped = numberLastSwapped;
1921  result.modifyCosts = modifyCosts;
1922}
1923/*
1924  Chooses incoming
1925  Puts flipped ones in list
1926  If necessary will modify costs
1927*/
1928void AbcSimplexDual::dualColumn2()
1929{
1930  CoinPartitionedVector &candidateList = usefulArray_[arrayForDualColumn_];
1931  CoinPartitionedVector &flipList = usefulArray_[arrayForFlipBounds_];
1932  //usefulArray_[arrayForTableauRow_].compact();
1933  usefulArray_[arrayForTableauRow_].computeNumberElements();
1934  if (candidateList.getNumElements() < 100)
1935    candidateList.compact();
1936  //usefulArray_[arrayForTableauRow_].compact();
1937  // But factorizations complain if <1.0e-8
1938  //acceptablePivot=CoinMax(acceptablePivot,1.0e-8);
1939  // do ratio test for normal iteration
1940  currentAcceptablePivot_ = 1.0e-1 * acceptablePivot_;
1941  if (numberIterations_ > 100)
1942    currentAcceptablePivot_ = acceptablePivot_;
1943  if (abcFactorization_->pivots() > 10 || (abcFactorization_->pivots() && sumDualInfeasibilities_))
1944    currentAcceptablePivot_ = 1.0e+3 * acceptablePivot_; // if we have iterated be more strict
1945  else if (abcFactorization_->pivots() > 5)
1946    currentAcceptablePivot_ = 1.0e+2 * acceptablePivot_; // if we have iterated be slightly more strict
1947  else if (abcFactorization_->pivots())
1948    currentAcceptablePivot_ = acceptablePivot_; // relax
1949  // If very large infeasibility use larger acceptable pivot
1950  if (abcFactorization_->pivots()) {
1951    double value = CoinMax(currentAcceptablePivot_, 1.0e-12 * fabs(dualOut_));
1952    currentAcceptablePivot_ = CoinMin(1.0e2 * currentAcceptablePivot_, value);
1953  }
1954  // return at once if no free but there were free
1955  if (lastFirstFree_ >= 0 && sequenceIn_ < 0) {
1956    candidateList.clearAndReset();
1957    upperTheta_ = COIN_DBL_MAX;
1958  }
1959  if (upperTheta_ == COIN_DBL_MAX) {
1960    usefulArray_[arrayForTableauRow_].compact();
1961    assert(!candidateList.getNumElements());
1962    return;
1963  }
1964  int numberSwapped = 0;
1965  //int numberRemaining = candidateList.getNumElements();
1966
1967  double useThru = 0.0; // for when variables flip
1968
1969  // If we think we need to modify costs (not if something from broad sweep)
1970  bool modifyCosts = false;
1971  // pivot elements
1972  //double *  COIN_RESTRICT array=candidateList.denseVector();
1973  // indices
1974  const double *COIN_RESTRICT abcLower = abcLower_;
1975  const double *COIN_RESTRICT abcUpper = abcUpper_;
1976  const double *COIN_RESTRICT abcSolution = abcSolution_;
1977  const double *COIN_RESTRICT abcDj = abcDj_;
1978  //const double *  COIN_RESTRICT abcPerturbation = abcPerturbation_;
1979  const unsigned char *COIN_RESTRICT internalStatus = internalStatus_;
1980  // We can also see if infeasible or pivoting on free
1981  const double multiplier[] = { 1.0, -1.0 };
1982
1983  // For list of flipped
1984  int numberLastSwapped = 0;
1985  int *COIN_RESTRICT flippedIndex = flipList.getIndices();
1986  double *COIN_RESTRICT flippedElement = flipList.denseVector();
1987  // If sum of bad small pivots too much
1988#ifdef MORE_CAREFUL
1989  bool badSumPivots = false;
1990#endif
1991  if (sequenceIn_ < 0) {
1992    int lastSequence = -1;
1993    double lastPivotValue = 0.0;
1994    double thisPivotValue = 0.0;
1995    dualColumnResult result;
1996    dualColumn2Most(result);
1997    sequenceIn_ = result.sequence;
1998    lastSequence = result.lastSequence;
1999    thisPivotValue = result.thisPivotValue;
2000    lastPivotValue = result.lastPivotValue;
2001    numberSwapped = result.numberSwapped;
2002    numberLastSwapped = result.numberLastSwapped;
2003    modifyCosts = result.modifyCosts;
2004    useThru = result.useThru;
2005#ifdef PRINT_RATIO_PROGRESS
2006    printf("Iteration %d - out of loop - sequenceIn %d pivotValue %g\n",
2007      numberIterations_, sequenceIn_, thisPivotValue);
2008#endif
2009    // can get here without sequenceIn_ set but with lastSequence
2010    if (sequenceIn_ < 0 && lastSequence >= 0) {
2011      // back to previous one
2012      sequenceIn_ = lastSequence;
2013      thisPivotValue = lastPivotValue;
2014#ifdef PRINT_RATIO_PROGRESS
2015      printf("BUT back one - out of loop - sequenceIn %d pivotValue %g\n",
2016        sequenceIn_, thisPivotValue);
2017#endif
2018    }
2019
2020    // Movement should be minimum for anti-degeneracy - unless
2021    // fixed variable out
2022    // no need if just one candidate
2023    assert(numberSwapped - numberLastSwapped >= 0);
2024    double minimumTheta;
2025    if (upperOut_ > lowerOut_ && numberSwapped - numberLastSwapped > 1)
2026#ifndef TRY_ABC_GUS
2027      minimumTheta = MINIMUMTHETA;
2028#else
2029      minimumTheta = minimumThetaMovement_;
2030#endif
2031    else
2032      minimumTheta = 0.0;
2033    if (sequenceIn_ >= 0) {
2034      double oldValue;
2035      alpha_ = thisPivotValue;
2036      // correct sign
2037      int iStatus = internalStatus[sequenceIn_] & 3;
2038      alpha_ *= multiplier[iStatus];
2039      oldValue = abcDj[sequenceIn_];
2040      theta_ = CoinMax(oldValue / alpha_, 0.0);
2041      if (theta_ < minimumTheta && fabs(alpha_) < 1.0e5 && 1) {
2042        // can't pivot to zero
2043        theta_ = minimumTheta;
2044      }
2045      // may need to adjust costs so all dual feasible AND pivoted is exactly 0
2046      if (modifyCosts
2047#ifdef MORE_CAREFUL
2048        && !badSumPivots
2049#endif
2050      ) {
2051        for (int i = numberLastSwapped; i < numberSwapped; i++) {
2052          int iSequence = flippedIndex[i];
2053          int iStatus = internalStatus[iSequence] & 3;
2054          // treat as if at lower bound
2055          double mult = multiplier[iStatus];
2056          double alpha = flippedElement[i];
2057          flippedElement[i] = 0.0;
2058          if (iSequence == sequenceIn_)
2059            continue;
2060          double oldValue = abcDj[iSequence] * mult;
2061          double value = oldValue - theta_ * alpha;
2062          //double perturbedTolerance = abcPerturbation[iSequence];
2063          if (-value > currentDualTolerance_) {
2064            //thisIncrease = true;
2065#if MODIFYCOST
2066            if (numberIterations_ >= normalDualColumnIteration_) {
2067              // modify cost to hit new tolerance
2068              double modification = alpha * theta_ - oldValue - currentDualTolerance_;
2069              if ((specialOptions_ & (2048 + 4096)) != 0) {
2070                if ((specialOptions_ & 2048) != 0) {
2071                  if (fabs(modification) < 1.0e-10)
2072                    modification = 0.0;
2073                } else {
2074                  if (fabs(modification) < 1.0e-12)
2075                    modification = 0.0;
2076                }
2077              }
2078              abcDj_[iSequence] += mult * modification;
2079              abcCost_[iSequence] += mult * modification;
2080              if (modification)
2081                numberChanged_++; // Say changed costs
2082            }
2083#endif
2084          }
2085        }
2086        numberSwapped = numberLastSwapped;
2087      }
2088    }
2089  }
2090
2091#ifdef MORE_CAREFUL
2092  // If we have done pivots and things look bad set alpha_ 0.0 to force factorization
2093  if ((badSumPivots || fabs(theta_ * badFree) > 10.0 * currentDualTolerance_) && abcFactorization_->pivots()) {
2094#if ABC_NORMAL_DEBUG > 0
2095    if (handler_->logLevel() > 1)
2096      printf("forcing re-factorization\n");
2097#endif
2098    sequenceIn_ = -1;
2099  }
2100#endif
2101  if (sequenceIn_ >= 0) {
2102    lowerIn_ = abcLower[sequenceIn_];
2103    upperIn_ = abcUpper[sequenceIn_];
2104    valueIn_ = abcSolution[sequenceIn_];
2105    dualIn_ = abcDj_[sequenceIn_];
2106
2107    if (numberTimesOptimal_) {
2108      // can we adjust cost back closer to original
2109      //*** add coding
2110    }
2111    if (numberIterations_ >= normalDualColumnIteration_) {
2112#if MODIFYCOST > 1
2113      // modify cost to hit zero exactly
2114      // so (dualIn_+modification)==theta_*alpha_
2115      //double perturbedTolerance = abcPerturbation[sequenceIn_];
2116      double modification = theta_ * alpha_ - (dualIn_ + goToTolerance2());
2117      // But should not move objective too much ??
2118#ifdef DONT_MOVE_OBJECTIVE
2119      double moveObjective = fabs(modification * abcSolution[sequenceIn_]);
2120      double smallMove = CoinMax(fabs(rawObjectiveValue_), 1.0e-3);
2121      if (moveObjective > smallMove) {
2122#if ABC_NORMAL_DEBUG > 0
2123        if (handler_->logLevel() > 1)
2124          printf("would move objective by %g - original mod %g sol value %g\n", moveObjective,
2125            modification, abcSolution[sequenceIn_]);
2126#endif
2127        modification *= smallMove / moveObjective;
2128      }
2129#endif
2130#ifdef MORE_CAREFUL
2131      if (badSumPivots)
2132        modification = 0.0;
2133#endif
2134      if (atFakeBound(sequenceIn_) && fabs(modification) < currentDualTolerance_)
2135        modification = 0.0;
2136      if ((specialOptions_ & (2048 + 4096)) != 0) {
2137        if ((specialOptions_ & 16384) != 0) {
2138          // in fast dual
2139          if (fabs(modification) < 1.0e-7)
2140            modification = 0.0;
2141        } else if ((specialOptions_ & 2048) != 0) {
2142          if (fabs(modification) < 1.0e-10)
2143            modification = 0.0;
2144        } else {
2145          if (fabs(modification) < 1.0e-12)
2146            modification = 0.0;
2147        }
2148      }
2149
2150      dualIn_ += modification;
2151      abcDj_[sequenceIn_] = dualIn_;
2152      abcCost_[sequenceIn_] += modification;
2153      if (modification)
2154        numberChanged_++; // Say changed costs
2155        //int costOffset = numberRows_+numberColumns_;
2156        //abcCost[sequenceIn_+costOffset] += modification; // save change
2157        //assert (fabs(modification)<1.0e-6);
2158#if ABC_NORMAL_DEBUG > 3
2159      if ((handler_->logLevel() & 32) && fabs(modification) > 1.0e-15)
2160        printf("exact %d new cost %g, change %g\n", sequenceIn_,
2161          abcCost_[sequenceIn_], modification);
2162#endif
2163#endif
2164    }
2165
2166    if (alpha_ < 0.0) {
2167      // as if from upper bound
2168      directionIn_ = -1;
2169      upperIn_ = valueIn_;
2170    } else {
2171      // as if from lower bound
2172      directionIn_ = 1;
2173      lowerIn_ = valueIn_;
2174    }
2175  } else {
2176    // no pivot
2177    alpha_ = 0.0;
2178  }
2179  // clear array
2180  //printf("dualout %g theta %g alpha %g dj %g thru %g swapped %d product %g guess %g useThru %g\n",
2181  //     dualOut_,theta_,alpha_,abcDj_[sequenceIn_],useThru,numberSwapped,
2182  //     fabs(dualOut_)*theta_,objectiveValue()+(fabs(dualOut_)-useThru)*theta_,useThru);
2183  objectiveChange_ = (fabs(dualOut_) - useThru) * theta_;
2184  rawObjectiveValue_ += objectiveChange_;
2185  //CoinZeroN(array, maximumNumberInList);
2186  //candidateList.setNumElements(0);
2187  candidateList.clearAndReset();
2188  flipList.setNumElements(numberSwapped);
2189  if (numberSwapped)
2190    flipList.setPackedMode(true);
2191#ifdef PRINT_RATIO_PROGRESS
2192  printf("%d flipped\n", numberSwapped);
2193#endif
2194}
2195/* Checks if tentative optimal actually means unbounded
2196   Returns -3 if not, 2 if is unbounded */
2197int AbcSimplexDual::checkUnbounded(CoinIndexedVector &ray,
2198  double changeCost)
2199{
2200  int status = 2; // say unbounded
2201  abcFactorization_->updateColumn(ray);
2202  // get reduced cost
2203  int number = ray.getNumElements();
2204  int *COIN_RESTRICT index = ray.getIndices();
2205  double *COIN_RESTRICT array = ray.denseVector();
2206  const int *COIN_RESTRICT abcPivotVariable = abcPivotVariable_;
2207  for (int i = 0; i < number; i++) {
2208    int iRow = index[i];
2209    int iPivot = abcPivotVariable[iRow];
2210    changeCost -= cost(iPivot) * array[iRow];
2211  }
2212  double way;
2213  if (changeCost > 0.0) {
2214    //try going down
2215    way = 1.0;
2216  } else if (changeCost < 0.0) {
2217    //try going up
2218    way = -1.0;
2219  } else {
2220#if ABC_NORMAL_DEBUG > 3
2221    printf("can't decide on up or down\n");
2222#endif
2223    way = 0.0;
2224    status = -3;
2225  }
2226  double movement = 1.0e10 * way; // some largish number
2227  double zeroTolerance = 1.0e-14 * dualBound_;
2228  for (int i = 0; i < number; i++) {
2229    int iRow = index[i];
2230    int iPivot = abcPivotVariable[iRow];
2231    double arrayValue = array[iRow];
2232    if (fabs(arrayValue) < zeroTolerance)
2233      arrayValue = 0.0;
2234    double newValue = solution(iPivot) + movement * arrayValue;
2235    if (newValue > upper(iPivot) + primalTolerance_ || newValue < lower(iPivot) - primalTolerance_)
2236      status = -3; // not unbounded
2237  }
2238  if (status == 2) {
2239    // create ray
2240    delete[] ray_;
2241    ray_ = new double[numberColumns_];
2242    CoinZeroN(ray_, numberColumns_);
2243    for (int i = 0; i < number; i++) {
2244      int iRow = index[i];
2245      int iPivot = abcPivotVariable[iRow];
2246      double arrayValue = array[iRow];
2247      if (iPivot < numberColumns_ && fabs(arrayValue) >= zeroTolerance)
2248        ray_[iPivot] = way * array[iRow];
2249    }
2250  }
2251  ray.clear();
2252  return status;
2253}
2254// Saves good status etc
2255void AbcSimplex::saveGoodStatus()
2256{
2257  // save arrays
2258  CoinMemcpyN(internalStatus_, numberTotal_, internalStatusSaved_);
2259  // save costs
2260  CoinMemcpyN(abcCost_, numberTotal_, abcCost_ + 2 * numberTotal_);
2261  CoinMemcpyN(abcSolution_, numberTotal_, solutionSaved_);
2262}
2263// Restores previous good status and says trouble
2264void AbcSimplex::restoreGoodStatus(int type)
2265{
2266#if PARTITION_ROW_COPY == 1
2267  stateOfProblem_ |= NEED_BASIS_SORT;
2268#endif
2269  if (type) {
2270    if (alphaAccuracy_ != -1.0)
2271      alphaAccuracy_ = -2.0;
2272    // trouble - restore solution
2273    CoinMemcpyN(internalStatusSaved_, numberTotal_, internalStatus_);
2274    CoinMemcpyN(abcCost_ + 2 * numberTotal_, numberTotal_, abcCost_);
2275    CoinMemcpyN(solutionSaved_, numberTotal_, abcSolution_);
2276    int *COIN_RESTRICT abcPivotVariable = abcPivotVariable_;
2277    // redo pivotVariable
2278    int numberBasic = 0;
2279    for (int i = 0; i < numberTotal_; i++) {
2280      if (getInternalStatus(i) == basic)
2281        abcPivotVariable[numberBasic++] = i;
2282    }
2283    assert(numberBasic == numberRows_);
2284    forceFactorization_ = 1; // a bit drastic but ..
2285    changeMade_++; // say something changed
2286  }
2287}
2288static int computePrimalsAndCheck(AbcSimplexDual *model, const int whichArray[2])
2289{
2290  CoinIndexedVector *array1 = model->usefulArray(whichArray[0]);
2291  CoinIndexedVector *array2 = model->usefulArray(whichArray[1]);
2292  int numberRefinements = model->computePrimals(array1, array2);
2293  model->checkPrimalSolution(true);
2294  return numberRefinements;
2295}
2296static int computeDualsAndCheck(AbcSimplexDual *model, const int whichArray[2])
2297{
2298  CoinIndexedVector *array1 = model->usefulArray(whichArray[0]);
2299  CoinIndexedVector *array2 = model->usefulArray(whichArray[1]);
2300  int numberRefinements = model->computeDuals(NULL, array1, array2);
2301  model->checkDualSolutionPlusFake();
2302  return numberRefinements;
2303}
2304// Computes solutions - 1 do duals, 2 do primals, 3 both
2305int AbcSimplex::gutsOfSolution(int type)
2306{
2307  double oldLargestPrimalError = largestPrimalError_;
2308  double oldLargestDualError = largestDualError_;
2309  AbcSimplexDual *dual = reinterpret_cast< AbcSimplexDual * >(this);
2310  // do work and check
2311  int numberRefinements = 0;
2312#if ABC_PARALLEL
2313  if (type != 3) {
2314#endif
2315    if ((type & 2) != 0) {
2316      //work space
2317      int whichArray[2];
2318      for (int i = 0; i < 2; i++)
2319        whichArray[i] = getAvailableArray();
2320      numberRefinements = computePrimalsAndCheck(dual, whichArray);
2321      for (int i = 0; i < 2; i++)
2322        setAvailableArray(whichArray[i]);
2323    }
2324    if ((type & 1) != 0
2325#if CAN_HAVE_ZERO_OBJ > 1
2326      && (specialOptions_ & 2097152) == 0
2327#endif
2328    ) {
2329      //work space
2330      int whichArray[2];
2331      for (int i = 0; i < 2; i++)
2332        whichArray[i] = getAvailableArray();
2333      numberRefinements += computeDualsAndCheck(dual, whichArray);
2334      for (int i = 0; i < 2; i++)
2335        setAvailableArray(whichArray[i]);
2336    }
2337#if ABC_PARALLEL
2338  } else {
2339#ifndef NDEBUG
2340    abcFactorization_->checkMarkArrays();
2341#endif
2342    // redo to do in parallel
2343    //work space
2344    int whichArray[5];
2345    for (int i = 1; i < 5; i++)
2346      whichArray[i] = getAvailableArray();
2347    if (parallelMode_ == 0) {
2348      numberRefinements += computePrimalsAndCheck(dual, whichArray + 3);
2349      numberRefinements += computeDualsAndCheck(dual, whichArray + 1);
2350    } else {
2351#if ABC_PARALLEL == 1
2352      threadInfo_[0].stuff[0] = whichArray[1];
2353      threadInfo_[0].stuff[1] = whichArray[2];
2354      stopStart_ = 1 + 32 * 1;
2355      startParallelStuff(1);
2356#else
2357      whichArray[0] = 1;
2358      CoinThreadInfo info;
2359      info.status = 1;
2360      info.stuff[0] = whichArray[1];
2361      info.stuff[1] = whichArray[2];
2362      int n = cilk_spawn computeDualsAndCheck(dual, whichArray + 1);
2363#endif
2364      numberRefinements = computePrimalsAndCheck(dual, whichArray + 3);
2365#if ABC_PARALLEL == 1
2366      numberRefinements += stopParallelStuff(1);
2367#else
2368      cilk_sync;
2369      numberRefinements += n;
2370#endif
2371    }
2372    for (int i = 1; i < 5; i++)
2373      setAvailableArray(whichArray[i]);
2374  }
2375#endif
2376  rawObjectiveValue_ += sumNonBasicCosts_;
2377  setClpSimplexObjectiveValue();
2378  if (handler_->logLevel() > 3 || (largestPrimalError_ > 1.0e-2 || largestDualError_ > 1.0e-2))
2379    handler_->message(CLP_SIMPLEX_ACCURACY, messages_)
2380      << largestPrimalError_
2381      << largestDualError_
2382      << CoinMessageEol;
2383  if ((largestPrimalError_ > 1.0e1 || largestDualError_ > 1.0e1)
2384    && numberRows_ > 100 && numberIterations_) {
2385#if ABC_NORMAL_DEBUG > 0
2386    printf("Large errors - primal %g dual %g\n",
2387      largestPrimalError_, largestDualError_);
2388#endif
2389    // Change factorization tolerance
2390    //if (abcFactorization_->zeroTolerance() > 1.0e-18)
2391    //abcFactorization_->zeroTolerance(1.0e-18);
2392  }
2393  bool notChanged = true;
2394  if (numberIterations_ && (forceFactorization_ > 2 || forceFactorization_ < 0 || abcFactorization_->pivotTolerance() < 0.9899999999) && (oldLargestDualError || oldLargestPrimalError)) {
2395    double useOldDualError = oldLargestDualError;
2396    double useDualError = largestDualError_;
2397    if (algorithm_ > 0 && abcNonLinearCost_ && abcNonLinearCost_->sumInfeasibilities()) {
2398      double factor = CoinMax(1.0, CoinMin(1.0e3, infeasibilityCost_ * 1.0e-6));
2399      useOldDualError /= factor;
2400      useDualError /= factor;
2401    }
2402    if ((largestPrimalError_ > 1.0e3 && oldLargestPrimalError * 1.0e2 < largestPrimalError_) || (useDualError > 1.0e3 && useOldDualError * 1.0e2 < useDualError)) {
2403      double pivotTolerance = abcFactorization_->pivotTolerance();
2404      double factor = (largestPrimalError_ > 1.0e10 || largestDualError_ > 1.0e10)
2405        ? 2.0
2406        : 1.2;
2407      if (pivotTolerance < 0.1)
2408        abcFactorization_->pivotTolerance(0.1);
2409      else if (pivotTolerance < 0.98999999)
2410        abcFactorization_->pivotTolerance(CoinMin(0.99, pivotTolerance * factor));
2411      notChanged = pivotTolerance == abcFactorization_->pivotTolerance();
2412#if defined(CLP_USEFUL_PRINTOUT) && !defined(GCC_4_9)
2413      if (pivotTolerance < 0.9899999) {
2414        double newTolerance = abcFactorization_->pivotTolerance();
2415        printf("Changing pivot tolerance from %g to %g and backtracking\n",
2416          pivotTolerance, newTolerance);
2417      }
2418      printf("because old,new primal error %g,%g - dual %g,%g pivot_tol %g\n",
2419        oldLargestPrimalError, largestPrimalError_,
2420        oldLargestDualError, largestDualError_,
2421        pivotTolerance);
2422#endif
2423      if (pivotTolerance < 0.9899999) {
2424        //largestPrimalError_=0.0;
2425        //largestDualError_=0.0;
2426        //returnCode=1;
2427      }
2428    }
2429  }
2430  if (progress_.iterationNumber_[0] > 0 && progress_.iterationNumber_[CLP_PROGRESS - 1] - progress_.iterationNumber_[0] < CLP_PROGRESS * 3 && abcFactorization_->pivotTolerance() < 0.25 && notChanged) {
2431    double pivotTolerance = abcFactorization_->pivotTolerance();
2432    abcFactorization_->pivotTolerance(pivotTolerance * 1.5);
2433#if defined(CLP_USEFUL_PRINTOUT) && !defined(GCC_4_9)
2434    double newTolerance = abcFactorization_->pivotTolerance();
2435    printf("Changing pivot tolerance from %g to %g - inverting too often\n",
2436      pivotTolerance, newTolerance);
2437#endif
2438  }
2439  return numberRefinements;
2440}
2441// Make non free variables dual feasible by moving to a bound
2442int AbcSimplexDual::makeNonFreeVariablesDualFeasible(bool changeCosts)
2443{
2444  const double multiplier[] = { 1.0, -1.0 }; // sign?????????????
2445  double dualT = -100.0 * currentDualTolerance_;
2446  int numberMoved = 0;
2447  int numberChanged = 0;
2448  double *COIN_RESTRICT solution = abcSolution_;
2449  double *COIN_RESTRICT lower = abcLower_;
2450  double *COIN_RESTRICT upper = abcUpper_;
2451  double saveNonBasicCosts = sumNonBasicCosts_;
2452  double largestCost = dualTolerance_;
2453  for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
2454    unsigned char iStatus = internalStatus_[iSequence] & 7;
2455    if (iStatus < 4) {
2456      double mult = multiplier[iStatus];
2457      double dj = abcDj_[iSequence] * mult;
2458      if (dj < dualT) {
2459        largestCost = CoinMax(fabs(abcPerturbation_[iSequence]), largestCost);
2460      }
2461    } else if (iStatus != 7) {
2462      largestCost = CoinMax(fabs(abcPerturbation_[iSequence]), largestCost);
2463    }
2464  }
2465  for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
2466    unsigned char iStatus = internalStatus_[iSequence] & 7;
2467    if (iStatus < 4) {
2468      double mult = multiplier[iStatus];
2469      double dj = abcDj_[iSequence] * mult;
2470      // ? should use perturbationArray
2471      if (dj < dualT && changeCosts) {
2472        double gap = upper[iSequence] - lower[iSequence];
2473        if (gap > 1.1 * currentDualBound_) {
2474          assert(getFakeBound(iSequence) == AbcSimplexDual::noFake);
2475          if (!iStatus) {
2476            upper[iSequence] = lower[iSequence] + currentDualBound_;
2477            setFakeBound(iSequence, AbcSimplexDual::upperFake);
2478          } else {
2479            lower[iSequence] = upper[iSequence] - currentDualBound_;
2480            setFakeBound(iSequence, AbcSimplexDual::lowerFake);
2481          }
2482        }
2483        double costDifference = abcPerturbation_[iSequence] - abcCost_[iSequence];
2484        if (costDifference * mult > 0.0) {
2485#if ABC_NORMAL_DEBUG > 0
2486          //printf("can improve situation on %d by %g\n",iSequence,costDifference*mult);
2487#endif
2488          dj += costDifference * mult;
2489          abcDj_[iSequence] += costDifference;
2490          abcCost_[iSequence] = abcPerturbation_[iSequence];
2491          costDifference = 0.0;
2492        }
2493        double changeInObjective = -gap * dj;
2494        // probably need to be more intelligent
2495        if ((changeInObjective > 1.0e3 || gap > CoinMax(0.5 * dualBound_, 1.0e3)) && dj > -1.0e-1) {
2496          // change cost
2497          costDifference = (-dj - 0.5 * dualTolerance_) * mult;
2498          if (fabs(costDifference) < 100000.0 * dualTolerance_) {
2499            dj += costDifference * mult;
2500            abcDj_[iSequence] += costDifference;
2501            abcCost_[iSequence] += costDifference;
2502            numberChanged++;
2503          }
2504        }
2505      }
2506      if (dj < dualT && !changeCosts) {
2507        numberMoved++;
2508        if (iStatus == 0) {
2509          // at lower bound
2510          // to upper bound
2511          setInternalStatus(iSequence, atUpperBound);
2512          solution[iSequence] = upper[iSequence];
2513          sumNonBasicCosts_ += abcCost_[iSequence] * (upper[iSequence] - lower[iSequence]);
2514        } else {
2515          // at upper bound
2516          // to lower bound
2517          setInternalStatus(iSequence, atLowerBound);
2518          solution[iSequence] = lower[iSequence];
2519          sumNonBasicCosts_ -= abcCost_[iSequence] * (upper[iSequence] - lower[iSequence]);
2520        }
2521      }
2522    }
2523  }
2524#if ABC_NORMAL_DEBUG > 0
2525  if (numberMoved + numberChanged)
2526    printf("makeDualfeasible %d moved, %d had costs changed - sum dual %g\n",
2527      numberMoved, numberChanged, sumDualInfeasibilities_);
2528#endif
2529  rawObjectiveValue_ += (sumNonBasicCosts_ - saveNonBasicCosts);
2530  setClpSimplexObjectiveValue();
2531  return numberMoved;
2532}
2533int AbcSimplexDual::whatNext()
2534{
2535  int returnCode = 0;
2536  assert(!stateOfIteration_);
2537  // see if update stable
2538#if ABC_NORMAL_DEBUG > 3
2539  if ((handler_->logLevel() & 32))
2540    printf("btran alpha %g, ftran alpha %g\n", btranAlpha_, alpha_);
2541#endif
2542  int numberPivots = abcFactorization_->pivots();
2543  //double checkValue = 1.0e-7; // numberPivots<5 ? 1.0e-7 : 1.0e-6;
2544  double checkValue = numberPivots ? 1.0e-7 : 1.0e-5;
2545  // relax if free variable in
2546  if ((internalStatus_[sequenceIn_] & 7) > 1)
2547    checkValue = 1.0e-5;
2548  // if can't trust much and long way from optimal then relax
2549  if (largestPrimalError_ > 10.0)
2550    checkValue = CoinMin(1.0e-4, 1.0e-8 * largestPrimalError_);
2551  if (fabs(btranAlpha_) < 1.0e-12 || fabs(alpha_) < 1.0e-12 || fabs(btranAlpha_ - alpha_) > checkValue * (1.0 + fabs(alpha_))) {
2552    handler_->message(CLP_DUAL_CHECK, messages_)
2553      << btranAlpha_
2554      << alpha_
2555      << CoinMessageEol;
2556    // see with more relaxed criterion
2557    double test;
2558    if (fabs(btranAlpha_) < 1.0e-8 || fabs(alpha_) < 1.0e-8)
2559      test = 1.0e-1 * fabs(alpha_);
2560    else
2561      test = 10.0 * checkValue; //1.0e-4 * (1.0 + fabs(alpha_));
2562    bool needFlag = (fabs(btranAlpha_) < 1.0e-12 || fabs(alpha_) < 1.0e-12 || fabs(btranAlpha_ - alpha_) > test);
2563#if ABC_NORMAL_DEBUG > 0
2564    if (fabs(btranAlpha_ + alpha_) < 1.0e-5 * (1.0 + fabs(alpha_))) {
2565      printf("alpha diff (%g,%g) %g check %g\n", btranAlpha_, alpha_, fabs(btranAlpha_ - alpha_), checkValue * (1.0 + fabs(alpha_)));
2566    }
2567#endif
2568    if (numberPivots) {
2569      if (needFlag && numberPivots < 10) {
2570        // need to reject something
2571        assert(sequenceOut_ >= 0);
2572        char x = isColumn(sequenceOut_) ? 'C' : 'R';
2573        handler_->message(CLP_SIMPLEX_FLAG, messages_)
2574          << x << sequenceWithin(sequenceOut_)
2575          << CoinMessageEol;
2576#if ABC_NORMAL_DEBUG > 0
2577        printf("flag %d as alpha  %g %g - after %d pivots\n", sequenceOut_,
2578          btranAlpha_, alpha_, numberPivots);
2579#endif
2580        // Make safer?
2581        abcFactorization_->saferTolerances(1.0, -1.03);
2582        setFlagged(sequenceOut_);
2583        // so can't happen immediately again
2584        sequenceOut_ = -1;
2585        lastBadIteration_ = numberIterations_; // say be more cautious
2586      }
2587      // Make safer?
2588      //if (numberPivots<5)
2589      //abcFactorization_->saferTolerances (-0.99, -1.01);
2590      problemStatus_ = -2; // factorize now
2591      returnCode = -2;
2592      stateOfIteration_ = 2;
2593    } else {
2594      if (needFlag) {
2595        // need to reject something
2596        assert(sequenceOut_ >= 0);
2597        char x = isColumn(sequenceOut_) ? 'C' : 'R';
2598        handler_->message(CLP_SIMPLEX_FLAG, messages_)
2599          << x << sequenceWithin(sequenceOut_)
2600          << CoinMessageEol;
2601#if ABC_NORMAL_DEBUG > 3
2602        printf("flag a %g %g\n", btranAlpha_, alpha_);
2603#endif
2604        // Make safer?
2605        double tolerance = abcFactorization_->pivotTolerance();
2606        abcFactorization_->saferTolerances(1.0, -1.03);
2607        setFlagged(sequenceOut_);
2608        // so can't happen immediately again
2609        sequenceOut_ = -1;
2610        //abcProgress_.clearBadTimes();
2611        lastBadIteration_ = numberIterations_; // say be more cautious
2612        if (fabs(alpha_) < 1.0e-10 && fabs(btranAlpha_) < 1.0e-8 && numberIterations_ > 100) {
2613          //printf("I think should declare infeasible\n");
2614          problemStatus_ = 1;
2615          returnCode = 1;
2616          stateOfIteration_ = 2;
2617        } else {
2618          stateOfIteration_ = 1;
2619          if (tolerance < abcFactorization_->pivotTolerance()) {
2620            stateOfIteration_ = 2;
2621            returnCode = -2;
2622          }
2623        }
2624      }
2625    }
2626  }
2627  if (!stateOfIteration_) {
2628    // check update
2629    returnCode = abcFactorization_->checkReplacePart2(pivotRow_, btranAlpha_, alpha_, ftAlpha_);
2630  }
2631  return returnCode;
2632}
2633/* Checks if finished. 
2634   type 0 initial
2635   1 normal
2636   2 initial but we have solved problem before
2637   3 can skip factorization
2638   Updates status */
2639void AbcSimplexDual::statusOfProblemInDual(int type)
2640{
2641  // If lots of iterations then adjust costs if large ones
2642  int numberPivots = abcFactorization_->pivots();
2643#if 0
2644  printf("statusOf %d pivots type %d pstatus %d forceFac %d dont %d\n",
2645         numberPivots,type,problemStatus_,forceFactorization_,dontFactorizePivots_);
2646#endif
2647  int saveType = type;
2648  double lastSumDualInfeasibilities = sumDualInfeasibilities_;
2649  double lastSumPrimalInfeasibilities = sumPrimalInfeasibilities_;
2650  if (type >= 4) {
2651    // force factorization
2652    type &= 3;
2653    numberPivots = 9999999;
2654  }
2655  //double realDualInfeasibilities = 0.0;
2656  double changeCost;
2657  bool unflagVariables = numberFlagged_ > 0;
2658  int weightsSaved = 0;
2659  int numberBlackMarks = 0;
2660  if (!type) {
2661    // be optimistic
2662    stateOfProblem_ &= ~PESSIMISTIC;
2663    // but use 0.1
2664    double newTolerance = CoinMax(0.1, saveData_.pivotTolerance_);
2665    abcFactorization_->pivotTolerance(newTolerance);
2666  }
2667  if (type != 3 && (numberPivots > dontFactorizePivots_ || numberIterations_ == baseIteration_)) {
2668    // factorize
2669    // save dual weights
2670    abcDualRowPivot_->saveWeights(this, 1);
2671    weightsSaved = 2;
2672    // is factorization okay?
2673    int factorizationStatus = 0;
2674    if ((largestPrimalError_ > 1.0e1 || largestDualError_ > 1.0e1)
2675      && numberRows_ > 100 && numberIterations_ > baseIteration_) {
2676      double newTolerance = CoinMin(1.1 * abcFactorization_->pivotTolerance(), 0.99);
2677      abcFactorization_->pivotTolerance(newTolerance);
2678    }
2679#ifdef EARLY_FACTORIZE
2680    if ((numberEarly_ & 0xffff) <= 5 || !usefulArray_[ABC_NUMBER_USEFUL - 1].getNumElements()) {
2681      //#define KEEP_TRYING_EARLY
2682#ifdef KEEP_TRYING_EARLY
2683      int savedEarly = numberEarly_ >> 16;
2684      int newEarly = (numberEarly_ & 0xffff);
2685      newEarly = CoinMin(savedEarly, 1 + 2 * newEarly);
2686      numberEarly_ = (savedEarly << 16) | newEarly;
2687#endif
2688      factorizationStatus = internalFactorize(type ? 1 : 2);
2689    } else {
2690      CoinIndexedVector &vector = usefulArray_[ABC_NUMBER_USEFUL - 1];
2691      bool bad = vector.getNumElements() < 0;
2692      const int *indices = vector.getIndices();
2693      double *array = vector.denseVector();
2694      int capacity = vector.capacity();
2695      capacity--;
2696      int numberPivotsStored = indices[capacity];
2697#if ABC_NORMAL_DEBUG > 0
2698      printf("early status %s nPivots %d\n", bad ? "bad" : "good", numberPivotsStored);
2699#endif
2700      if (!bad) {
2701        // do remaining pivots
2702        int iterationsDone = abcEarlyFactorization_->pivots();
2703        factorizationStatus = abcEarlyFactorization_->replaceColumns(this, vector,
2704          iterationsDone, numberPivotsStored, true);
2705        if (!factorizationStatus) {
2706          // should be able to compare e.g. basics 1.0
2707#define GOOD_EARLY
2708#ifdef GOOD_EARLY
2709          AbcSimplexFactorization *temp = abcFactorization_;
2710          abcFactorization_ = abcEarlyFactorization_;
2711          abcEarlyFactorization_ = temp;
2712#endif
2713          moveToBasic();
2714#ifndef GOOD_EARLY
2715          int *pivotSave = CoinCopyOfArray(abcPivotVariable_, numberRows_);
2716          factorizationStatus = internalFactorize(type ? 1 : 2);
2717          {
2718            double *array = usefulArray_[3].denseVector();
2719            double *array2 = usefulArray_[4].denseVector();
2720            for (int iPivot = 0; iPivot < numberRows_; iPivot++) {
2721              int iSequence = abcPivotVariable_[iPivot];
2722              unpack(usefulArray_[3], iSequence);
2723              unpack(usefulArray_[4], iSequence);
2724              abcFactorization_->updateColumn(usefulArray_[3]);
2725              abcEarlyFactorization_->updateColumn(usefulArray_[4]);
2726              assert(fabs(array[iPivot] - 1.0) < 1.0e-4);
2727              array[iPivot] = 0.0;
2728              int iPivot2;
2729              for (iPivot2 = 0; iPivot2 < numberRows_; iPivot2++) {
2730                if (pivotSave[iPivot2] == iSequence)
2731                  break;
2732              }
2733              assert(iPivot2 < numberRows_);
2734              assert(fabs(array2[iPivot2] - 1.0) < 1.0e-4);
2735              array2[iPivot2] = 0.0;
2736              for (int i = 0; i < numberRows_; i++)
2737                assert(fabs(array[i]) < 1.0e-4);
2738              for (int i = 0; i < numberRows_; i++)
2739                assert(fabs(array2[i]) < 1.0e-4);
2740              usefulArray_[3].clear();
2741              usefulArray_[4].clear();
2742            }
2743            delete[] pivotSave;
2744          }
2745#endif
2746        } else {
2747          bad = true;
2748        }
2749      }
2750      // clean up
2751      vector.setNumElements(0);
2752      for (int i = 0; i < 2 * numberPivotsStored; i++) {
2753        array[--capacity] = 0.0;
2754      }
2755      if (bad) {
2756        // switch off
2757#if ABC_NORMAL_DEBUG > 0
2758        printf("switching off early factorization(in statusOf)\n");
2759#endif
2760#ifndef KEEP_TRYING_EARLY
2761        numberEarly_ = 0;
2762        delete abcEarlyFactorization_;
2763        abcEarlyFactorization_ = NULL;
2764#else
2765        numberEarly_ &= 0xffff0000;
2766#endif
2767        factorizationStatus = internalFactorize(1);
2768      }
2769    }
2770#else
2771    factorizationStatus = internalFactorize(type ? 1 : 2);
2772#endif
2773    if (factorizationStatus) {
2774      if (type == 1 && lastFlaggedIteration_ != 1) {
2775        // use for real disaster
2776        lastFlaggedIteration_ = 1;
2777        // no - restore previous basis
2778        unflagVariables = false;
2779        // restore weights
2780        weightsSaved = 4;
2781        changeMade_++; // say something changed
2782        // Keep any flagged variables
2783        for (int i = 0; i < numberTotal_; i++) {
2784          if (flagged(i))
2785            internalStatusSaved_[i] |= 64; //say flagged
2786        }
2787        restoreGoodStatus(type);
2788        // get correct bounds on all variables
2789        resetFakeBounds(1);
2790        numberBlackMarks = 10;
2791        if (sequenceOut_ >= 0) {
2792          // need to reject something
2793          char x = isColumn(sequenceOut_) ? 'C' : 'R';
2794          handler_->message(CLP_SIMPLEX_FLAG, messages_)
2795            << x << sequenceWithin(sequenceOut_)
2796            << CoinMessageEol;
2797#if ABC_NORMAL_DEBUG > 3
2798          printf("flag d\n");
2799#endif
2800          setFlagged(sequenceOut_);
2801          lastBadIteration_ = numberIterations_; // say be more cautious
2802          // so can't happen immediately again
2803          sequenceOut_ = -1;
2804        }
2805        //abcProgress_.clearBadTimes();
2806
2807        // Go to safe
2808        abcFactorization_->pivotTolerance(0.99);
2809      } else {
2810#if ABC_NORMAL_DEBUG > 0
2811        printf("Bad initial basis\n");
2812#endif
2813      }
2814      if (internalFactorize(1)) {
2815        // try once more
2816        if (internalFactorize(0))
2817          abort();
2818      }
2819#if PARTITION_ROW_COPY == 1
2820      int iVector = getAvailableArray();
2821      abcMatrix_->sortUseful(usefulArray_[iVector]);
2822      setAvailableArray(iVector);
2823#endif
2824    } else {
2825      lastFlaggedIteration_ = 0;
2826    }
2827  }
2828  // get primal and dual solutions
2829  numberBlackMarks += gutsOfSolution(3) / 2;
2830  assert(type != 0 || numberIterations_ == baseIteration_); // may be wrong
2831  if (!type) {
2832    initialSumInfeasibilities_ = sumPrimalInfeasibilities_;
2833    initialNumberInfeasibilities_ = numberPrimalInfeasibilities_;
2834    CoinAbcMemcpy(solutionSaved_ + numberTotal_, abcSolution_, numberTotal_);
2835    CoinAbcMemcpy(internalStatusSaved_ + numberTotal_, internalStatus_, numberTotal_);
2836  }
2837  double savePrimalInfeasibilities = sumPrimalInfeasibilities_;
2838  if ((moreSpecialOptions_ & 16384) != 0 && (moreSpecialOptions_ & 8192) == 0) {
2839    if (numberIterations_ == baseIteration_) {
2840      // If primal feasible and looks suited to primal go to primal
2841      if (!sumPrimalInfeasibilities_ && sumDualInfeasibilities_ > 1.0 && numberRows_ * 4 > numberColumns_) {
2842        // do primal
2843        // mark so can perturb
2844        initialSumInfeasibilities_ = -1.0;
2845        problemStatus_ = 10;
2846        return;
2847      }
2848    } else if ((numberIterations_ > 1000 && numberIterations_ * 10 > numberRows_ && numberIterations_ * 4 < numberRows_) || (sumPrimalInfeasibilities_ + sumDualInfeasibilities_ > 1.0e10 * (initialSumInfeasibilities_ + 1000.0) && numberIterations_ > 10000 && sumDualInfeasibilities_ > 1.0)) {
2849      if (sumPrimalInfeasibilities_ + 1.0e2 * sumDualInfeasibilities_ > 2000.0 * (initialSumInfeasibilities_ + 100.0)) {
2850#if ABC_NORMAL_DEBUG > 0
2851        printf("sump %g sumd %g initial %g\n", sumPrimalInfeasibilities_, sumDualInfeasibilities_,
2852          initialSumInfeasibilities_);
2853#endif
2854        // but only if average gone up as well
2855        if (sumPrimalInfeasibilities_ * initialNumberInfeasibilities_ > 100.0 * (initialSumInfeasibilities_ + 10.0) * numberPrimalInfeasibilities_) {
2856#if ABC_NORMAL_DEBUG > 0
2857          printf("returning to do primal\n");
2858#endif
2859          // do primal
2860          problemStatus_ = 10;
2861          CoinAbcMemcpy(abcSolution_, solutionSaved_ + numberTotal_, numberTotal_);
2862          CoinAbcMemcpy(internalStatus_, internalStatusSaved_ + numberTotal_, numberTotal_);
2863          // mark so can perturb
2864          initialSumInfeasibilities_ = -1.0;
2865          return;
2866        }
2867      }
2868    }
2869  } else if ((sumPrimalInfeasibilities_ + sumDualInfeasibilities_ > 1.0e6 * (initialSumInfeasibilities_ + 1000.0) && (moreSpecialOptions_ & 8192) == 0 && numberIterations_ > 1000 && numberIterations_ * 10 > numberRows_ && numberIterations_ * 4 < numberRows_) || (sumPrimalInfeasibilities_ + sumDualInfeasibilities_ > 1.0e10 * (initialSumInfeasibilities_ + 1000.0) && (moreSpecialOptions_ & 8192) == 0 && numberIterations_ > 10000 && sumDualInfeasibilities_ > 1.0)) {
2870    // do primal
2871    problemStatus_ = 10;
2872    CoinAbcMemcpy(abcSolution_, solutionSaved_ + numberTotal_, numberTotal_);
2873    CoinAbcMemcpy(internalStatus_, internalStatusSaved_ + numberTotal_, numberTotal_);
2874    // mark so can perturb
2875    initialSumInfeasibilities_ = -1.0;
2876    return;
2877  }
2878  //double saveDualInfeasibilities=sumDualInfeasibilities_;
2879  // 0 don't ignore, 1 ignore errors, 2 values pass
2880  int ignoreStuff = (firstFree_ < 0 && lastFirstFree_ < 0) ? 0 : 1;
2881  if ((stateOfProblem_ & VALUES_PASS) != 0)
2882    ignoreStuff = 2;
2883  {
2884    double thisObj = rawObjectiveValue();
2885    double lastObj = abcProgress_.lastObjective(0);
2886    if ((lastObj > thisObj + 1.0e-3 * CoinMax(fabs(thisObj), fabs(lastObj)) + 1.0e-4 || (sumPrimalInfeasibilities_ > 10.0 * lastSumPrimalInfeasibilities + 1.0e3 && sumDualInfeasibilities_ > 10.0 * CoinMax(lastSumDualInfeasibilities, 1.0)))
2887      && !ignoreStuff) {
2888#if ABC_NORMAL_DEBUG > 0
2889      printf("bad - objective backwards %g %g errors %g %g suminf %g %g\n", lastObj, thisObj, largestDualError_, largestPrimalError_, sumDualInfeasibilities_, sumPrimalInfeasibilities_);
2890#endif
2891      numberBlackMarks++;
2892      if (largestDualError_ > 1.0e-3 || largestPrimalError_ > 1.0e-3 || (sumDualInfeasibilities_ > 10.0 * CoinMax(lastSumDualInfeasibilities, 1.0) && firstFree_ < 0)) {
2893#if ABC_NORMAL_DEBUG > 0
2894        printf("backtracking %d iterations - dual error %g primal error %g - sumDualInf %g\n",
2895          numberPivots,
2896          largestDualError_, largestPrimalError_, sumDualInfeasibilities_);
2897#endif
2898        // restore previous basis
2899        unflagVariables = false;
2900        // restore weights
2901        weightsSaved = 4;
2902        changeMade_++; // say something changed
2903        // Keep any flagged variables
2904        for (int i = 0; i < numberTotal_; i++) {
2905          if (flagged(i))
2906            internalStatusSaved_[i] |= 64; //say flagged
2907        }
2908        restoreGoodStatus(type);
2909        // get correct bounds on all variables
2910        resetFakeBounds(1);
2911        numberBlackMarks += 10;
2912        if (sequenceOut_ >= 0) {
2913          // need to reject something
2914#if ABC_NORMAL_DEBUG > 0
2915          if (flagged(sequenceOut_))
2916            printf("BAD %x already flagged\n", sequenceOut_);
2917#endif
2918          char x = isColumn(sequenceOut_) ? 'C' : 'R';
2919          handler_->message(CLP_SIMPLEX_FLAG, messages_)
2920            << x << sequenceWithin(sequenceOut_)
2921            << CoinMessageEol;
2922#if ABC_NORMAL_DEBUG > 3
2923          printf("flag d\n");
2924#endif
2925          setFlagged(sequenceOut_);
2926          lastBadIteration_ = numberIterations_; // say be more cautious
2927#if ABC_NORMAL_DEBUG > 0
2928        } else {
2929          printf("backtracking - no valid sequence out\n");
2930#endif
2931        }
2932        // so can't happen immediately again
2933        sequenceOut_ = -1;
2934        //abcProgress_.clearBadTimes();
2935        // Go to safe
2936        //abcFactorization_->pivotTolerance(0.99);
2937        if (internalFactorize(1)) {
2938          abort();
2939        }
2940        moveToBasic();
2941#if PARTITION_ROW_COPY == 1
2942        int iVector = getAvailableArray();
2943        abcMatrix_->sortUseful(usefulArray_[iVector]);
2944        setAvailableArray(iVector);
2945#endif
2946        // get primal and dual solutions
2947        numberBlackMarks += gutsOfSolution(3);
2948        //gutsOfSolution(1);
2949        makeNonFreeVariablesDualFeasible();
2950      }
2951    }
2952  }
2953  //#define PAN
2954#ifdef PAN
2955  if (!type && (specialOptions_ & COIN_CBC_USING_CLP) == 0 && firstFree_ < 0) {
2956    // first time - set large bad ones superbasic
2957    int whichArray = getAvailableArray();
2958    double *COIN_RESTRICT work = usefulArray_[whichArray].denseVector();
2959    int number = 0;
2960    int *COIN_RESTRICT which = usefulArray_[whichArray].getIndices();
2961    const double *COIN_RESTRICT reducedCost = abcDj_;
2962    const double *COIN_RESTRICT lower = abcLower_;
2963    const double *COIN_RESTRICT upper = abcUpper_;
2964#if PRINT_PAN
2965    int nAtUb = 0;
2966#endif
2967    for (int i = 0; i < numberTotal_; i++) {
2968      Status status = getInternalStatus(i);
2969      double djValue = reducedCost[i];
2970      double badValue = 0.0;
2971      if (status == atLowerBound && djValue < -currentDualTolerance_) {
2972        double gap = upper[i] - lower[i];
2973        if (gap > 100000000.0)
2974          badValue = -djValue * CoinMin(gap, 1.0e10);
2975      } else if (status == atUpperBound && djValue > currentDualTolerance_) {
2976        double gap = upper[i] - lower[i];
2977        if (gap > 100000000.0) {
2978          badValue = djValue * CoinMin(gap, 1.0e10);
2979#if PRINT_PAN
2980          nAtUb++;
2981#endif
2982        }
2983      }
2984      if (badValue) {
2985        work[number] = badValue;
2986        which[number++] = i;
2987      }
2988    }
2989    // for now don't do if primal feasible (probably should go to primal)
2990    if (!numberPrimalInfeasibilities_) {
2991      number = 0;
2992#if ABC_NORMAL_DEBUG > 1
2993      printf("XXXX doesn't work if primal feasible - also fix switch to Clp primal\n");
2994#endif
2995    }
2996    if (number && ignoreStuff != 2) {
2997      CoinSort_2(work, work + number, which);
2998      int numberToDo = CoinMin(number, numberRows_ / 10); // ??
2999#if PRINT_PAN > 1
3000      CoinSort_2(which, which + numberToDo, work);
3001      int n = 0;
3002#endif
3003#if PRINT_PAN
3004      printf("Panning %d variables (out of %d at lb and %d at ub)\n", numberToDo, number - nAtUb, nAtUb);
3005#endif
3006      firstFree_ = numberTotal_;
3007      for (int i = 0; i < numberToDo; i++) {
3008        int iSequence = which[i];
3009#if PRINT_PAN > 1
3010        if (!n)
3011          printf("Pan ");
3012        printf("%d ", iSequence);
3013        n++;
3014        if (n == 10) {
3015          n = 0;
3016          printf("\n");
3017        }
3018#endif
3019        setInternalStatus(iSequence, superBasic);
3020        firstFree_ = CoinMin(firstFree_, iSequence);
3021      }
3022#if PRINT_PAN > 1
3023      if (n)
3024        printf("\n");
3025#endif
3026      ordinaryVariables_ = 0;
3027      CoinAbcMemset0(work, number);
3028      stateOfProblem_ |= FAKE_SUPERBASIC;
3029    }
3030    setAvailableArray(whichArray);
3031  }
3032#endif
3033  if (!sumOfRelaxedPrimalInfeasibilities_ && sumDualInfeasibilities_ && type && saveType < 9
3034    && ignoreStuff != 2) {
3035    // say been feasible
3036    abcState_ |= CLP_ABC_BEEN_FEASIBLE;
3037    int numberFake = numberAtFakeBound();
3038    if (!numberFake) {
3039      makeNonFreeVariablesDualFeasible();
3040      numberDualInfeasibilities_ = 0;
3041      sumDualInfeasibilities_ = 0.0;
3042      // just primal
3043      gutsOfSolution(2);
3044    }
3045  }
3046  bool computeNewSumInfeasibilities = numberIterations_ == baseIteration_;
3047  if ((!type || (firstFree_ < 0 && lastFirstFree_ >= 0)) && ignoreStuff != 2) {
3048#ifdef PAN
3049    if ((stateOfProblem_ & FAKE_SUPERBASIC) != 0 && firstFree_ < 0 && lastFirstFree_ >= 0 && type) {
3050      // clean up
3051#if PRINT_PAN
3052      printf("Pan switch off after %d iterations\n", numberIterations_);
3053#endif
3054      computeNewSumInfeasibilities = true;
3055      ordinaryVariables_ = 1;
3056      const double *COIN_RESTRICT reducedCost = abcDj_;
3057      const double *COIN_RESTRICT lower = abcLower_;
3058      const double *COIN_RESTRICT upper = abcUpper_;
3059      for (int i = 0; i < numberTotal_; i++) {
3060        Status status = getInternalStatus(i);
3061        if (status == superBasic) {
3062          double djValue = reducedCost[i];
3063          double value = abcSolution_[i];
3064          if (value < lower[i] + primalTolerance_) {
3065            setInternalStatus(i, atLowerBound);
3066#if ABC_NORMAL_DEBUG > 0
3067            if (djValue < -dualTolerance_)
3068              printf("Odd at lb %d dj %g\n", i, djValue);
3069#endif
3070          } else if (value > upper[i] - primalTolerance_) {
3071            setInternalStatus(i, atUpperBound);
3072#if ABC_NORMAL_DEBUG > 0
3073            if (djValue > dualTolerance_)
3074              printf("Odd at ub %d dj %g\n", i, djValue);
3075#endif
3076          } else {
3077            assert(status != superBasic);
3078          }
3079        }
3080        if (status == isFree)
3081          ordinaryVariables_ = 0;
3082      }
3083      stateOfProblem_ &= ~FAKE_SUPERBASIC;
3084    }
3085#endif
3086    // set up perturbation and dualBound_
3087    // make dual feasible
3088    if (firstFree_ < 0) {
3089      if (numberFlagged_) {
3090        numberFlagged_ = 0;
3091        for (int iRow = 0; iRow < numberRows_; iRow++) {
3092          int iPivot = abcPivotVariable_[iRow];
3093          if (flagged(iPivot)) {
3094            clearFlagged(iPivot);
3095          }
3096        }
3097      }
3098      int numberChanged = resetFakeBounds(0);
3099#if ABC_NORMAL_DEBUG > 0
3100      if (numberChanged > 0)
3101        printf("Re-setting %d fake bounds %d gone dual infeasible - can expect backward objective\n",
3102          numberChanged, numberDualInfeasibilities_);
3103#endif
3104      if (numberDualInfeasibilities_ || numberChanged != -1) {
3105        makeNonFreeVariablesDualFeasible();
3106        numberDualInfeasibilities_ = 0;
3107        sumDualInfeasibilities_ = 0.0;
3108        sumOfRelaxedDualInfeasibilities_ = 0.0;
3109        // just primal
3110        gutsOfSolution(2);
3111        if (!numberIterations_ && sumPrimalInfeasibilities_ > 1.0e5 * (savePrimalInfeasibilities + 1.0e3) && (moreSpecialOptions_ & (256 | 8192)) == 0) {
3112          // Use primal
3113          problemStatus_ = 10;
3114          // mark so can perturb
3115          initialSumInfeasibilities_ = -1.0;
3116          CoinAbcMemcpy(abcSolution_, solutionSaved_, numberTotal_);
3117          CoinAbcMemcpy(internalStatus_, internalStatusSaved_, numberTotal_);
3118          return;
3119        }
3120        // redo
3121        if (computeNewSumInfeasibilities)
3122          initialSumInfeasibilities_ = sumPrimalInfeasibilities_;
3123        //computeObjective();
3124      }
3125    }
3126    //bounceTolerances(type);
3127  } else if (numberIterations_ > numberRows_ && stateDualColumn_ == -2) {
3128    bounceTolerances(4);
3129  }
3130  // If bad accuracy treat as singular
3131  if ((largestPrimalError_ > 1.0e15 || largestDualError_ > 1.0e15) && numberIterations_ > 1000) {
3132    problemStatus_ = 10; //abort();
3133#if ABC_NORMAL_DEBUG > 0
3134    printf("Big problems - errors %g %g try primal\n", largestPrimalError_, largestDualError_);
3135#endif
3136  } else if (goodAccuracy() && numberIterations_ > lastBadIteration_ + 200) {
3137    // Can reduce tolerance
3138    double newTolerance = CoinMax(0.995 * abcFactorization_->pivotTolerance(), saveData_.pivotTolerance_);
3139    if (!type)
3140      newTolerance = saveData_.pivotTolerance_;
3141    if (newTolerance > abcFactorization_->minimumPivotTolerance())
3142      abcFactorization_->pivotTolerance(newTolerance);
3143  }
3144  bestObjectiveValue_ = CoinMax(bestObjectiveValue_,
3145    rawObjectiveValue_ - bestPossibleImprovement_);
3146  // Double check infeasibility if no action
3147  if (abcProgress_.lastIterationNumber(0) == numberIterations_) {
3148    if (abcDualRowPivot_->looksOptimal()) {
3149      // say been feasible
3150      abcState_ |= CLP_ABC_BEEN_FEASIBLE;
3151      numberPrimalInfeasibilities_ = 0;
3152      sumPrimalInfeasibilities_ = 0.0;
3153    }
3154  } else {
3155    double thisObj = rawObjectiveValue_ - bestPossibleImprovement_;
3156#ifdef CLP_INVESTIGATE
3157    assert(bestPossibleImprovement_ > -1000.0 && rawObjectiveValue_ > -1.0e100);
3158    if (bestPossibleImprovement_)
3159      printf("obj %g add in %g -> %g\n", rawObjectiveValue_, bestPossibleImprovement_,
3160        thisObj);
3161#endif
3162    double lastObj = abcProgress_.lastObjective(0);
3163#ifndef NDEBUG
3164#if ABC_NORMAL_DEBUG > 3
3165    resetFakeBounds(-1);
3166#endif
3167#endif
3168    if (lastObj < thisObj - 1.0e-5 * CoinMax(fabs(thisObj), fabs(lastObj)) - 1.0e-3 && saveType < 9) {
3169      numberTimesOptimal_ = 0;
3170    }
3171  }
3172  // Up tolerance if looks a bit odd
3173  if (numberIterations_ > CoinMax(1000, numberRows_ >> 4) && (specialOptions_ & 64) != 0) {
3174    if (sumPrimalInfeasibilities_ && sumPrimalInfeasibilities_ < 1.0e5) {
3175      int backIteration = abcProgress_.lastIterationNumber(CLP_PROGRESS - 1);
3176      if (backIteration > 0 && numberIterations_ - backIteration < 9 * CLP_PROGRESS) {
3177        if (abcFactorization_->pivotTolerance() < 0.9) {
3178          // up tolerance
3179          abcFactorization_->pivotTolerance(CoinMin(abcFactorization_->pivotTolerance() * 1.05 + 0.02, 0.91));
3180          //printf("tol now %g\n",abcFactorization_->pivotTolerance());
3181          abcProgress_.clearIterationNumbers();
3182        }
3183      }
3184    }
3185  }
3186  // Check if looping
3187  int loop;
3188  if (type != 2)
3189    loop = abcProgress_.looping();
3190  else
3191    loop = -1;
3192  if (abcProgress_.reallyBadTimes() > 10) {
3193    problemStatus_ = 10; // instead - try other algorithm
3194#if ABC_NORMAL_DEBUG > 3
3195    printf("returning at %d\n", __LINE__);
3196#endif
3197  }
3198  if (loop >= 0) {
3199    problemStatus_ = loop; //exit if in loop
3200    if (!problemStatus_) {
3201      // declaring victory
3202      // say been feasible
3203      abcState_ |= CLP_ABC_BEEN_FEASIBLE;
3204      numberPrimalInfeasibilities_ = 0;
3205      sumPrimalInfeasibilities_ = 0.0;
3206    } else {
3207      problemStatus_ = 10; // instead - try other algorithm
3208#if ABC_NORMAL_DEBUG > 3
3209      printf("returning at %d\n", __LINE__);
3210#endif
3211    }
3212    return;
3213  }
3214  // really for free variables in
3215  if ((progressFlag_ & 2) != 0) {
3216    //situationChanged = 2;
3217  }
3218  progressFlag_ &= (~3); //reset progress flag
3219  // mark as having gone optimal if looks like it
3220  if (!numberPrimalInfeasibilities_ && !numberDualInfeasibilities_)
3221    progressFlag_ |= 8;
3222  if (handler_->detail(CLP_SIMPLEX_STATUS, messages_) < 100) {
3223    handler_->message(CLP_SIMPLEX_STATUS, messages_)
3224      << numberIterations_ << objectiveValue();
3225    handler_->printing(sumPrimalInfeasibilities_ > 0.0)
3226      << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
3227    handler_->printing(sumDualInfeasibilities_ > 0.0)
3228      << sumDualInfeasibilities_ << numberDualInfeasibilities_;
3229    handler_->printing(numberDualInfeasibilitiesWithoutFree_
3230      < numberDualInfeasibilities_)
3231      << numberDualInfeasibilitiesWithoutFree_;
3232    handler_->message() << CoinMessageEol;
3233  }
3234  double approximateObjective = rawObjectiveValue_;
3235  //realDualInfeasibilities = sumDualInfeasibilities_;
3236  // If we need to carry on cleaning variables
3237  if (!numberPrimalInfeasibilities_ && (specialOptions_ & 1024) != 0 && CLEAN_FIXED) {
3238    for (int iRow = 0; iRow < numberRows_; iRow++) {
3239      int iPivot = abcPivotVariable_[iRow];
3240      if (!flagged(iPivot) && pivoted(iPivot)) {
3241        // carry on
3242        numberPrimalInfeasibilities_ = -1;
3243        sumOfRelaxedPrimalInfeasibilities_ = 1.0;
3244        sumPrimalInfeasibilities_ = 1.0;
3245        break;
3246      }
3247    }
3248  }
3249  /* If we are primal feasible and any dual infeasibilities are on
3250     free variables then it is better to go to primal */
3251  if (!numberPrimalInfeasibilities_ && !numberDualInfeasibilitiesWithoutFree_ && numberDualInfeasibilities_) {
3252    // say been feasible
3253    abcState_ |= CLP_ABC_BEEN_FEASIBLE;
3254    if (type && ignoreStuff != 2)
3255      problemStatus_ = 10;
3256    else
3257      numberPrimalInfeasibilities_ = 1;
3258  }
3259  // check optimal
3260  // give code benefit of doubt
3261  if (sumOfRelaxedDualInfeasibilities_ == 0.0 && sumOfRelaxedPrimalInfeasibilities_ == 0.0 && ignoreStuff != 2) {
3262    // say been feasible
3263    abcState_ |= CLP_ABC_BEEN_FEASIBLE;
3264    // say optimal (with these bounds etc)
3265    numberDualInfeasibilities_ = 0;
3266    sumDualInfeasibilities_ = 0.0;
3267    numberPrimalInfeasibilities_ = 0;
3268    sumPrimalInfeasibilities_ = 0.0;
3269  } else if (sumOfRelaxedDualInfeasibilities_ > 1.0e-3 && firstFree_ < 0 && ignoreStuff != 2
3270    && numberIterations_ > lastCleaned_ + 10) {
3271#if ABC_NORMAL_DEBUG > 0
3272    printf("Making dual feasible\n");
3273#endif
3274    makeNonFreeVariablesDualFeasible(true);
3275    lastCleaned_ = numberIterations_;
3276    numberDualInfeasibilities_ = 0;
3277    sumDualInfeasibilities_ = 0.0;
3278    // just primal
3279    gutsOfSolution(2);
3280  }
3281  // see if cutoff reached
3282  double limit = 0.0;
3283  getDblParam(ClpDualObjectiveLimit, limit);
3284  int numberChangedBounds = 0;
3285  if (primalFeasible() && ignoreStuff != 2) {
3286    // may be optimal - or may be bounds are wrong
3287    numberChangedBounds = changeBounds(0, changeCost);
3288    if (numberChangedBounds)
3289      gutsOfSolution(2);
3290  }
3291  if (primalFeasible() && ignoreStuff != 2) {
3292    if (numberChangedBounds <= 0 && !numberDualInfeasibilities_) {
3293      //looks optimal - do we need to reset tolerance
3294      if (lastCleaned_ < numberIterations_ && /*numberTimesOptimal_*/ stateDualColumn_ < 4 && (numberChanged_ || (specialOptions_ & 4096) == 0)) {
3295#if CLP_CAN_HAVE_ZERO_OBJ
3296        if ((specialOptions_ & 2097152) == 0) {
3297#endif
3298          numberTimesOptimal_++;
3299          if (stateDualColumn_ == -2) {
3300#if ABC_NORMAL_DEBUG > 0
3301            printf("Going straight from state -2 to 0\n");
3302#endif
3303            stateDualColumn_ = -1;
3304          }
3305          changeMade_++; // say something changed
3306          if (numberTimesOptimal_ == 1) {
3307            if (fabs(currentDualTolerance_ - dualTolerance_) > 1.0e-12) {
3308#if ABC_NORMAL_DEBUG > 0
3309              printf("changing dual tolerance from %g to %g (numberTimesOptimal_==%d)\n",
3310                currentDualTolerance_, dualTolerance_, numberTimesOptimal_);
3311#endif
3312              currentDualTolerance_ = dualTolerance_;
3313            }
3314          } else {
3315            if (numberTimesOptimal_ == 2) {
3316              // better to have small tolerance even if slower
3317              abcFactorization_->zeroTolerance(CoinMin(abcFactorization_->zeroTolerance(), 1.0e-15));
3318            }
3319#if ABC_NORMAL_DEBUG > 0
3320            printf("changing dual tolerance from %g to %g (numberTimesOptimal_==%d)\n",
3321              currentDualTolerance_, dualTolerance_ * pow(2.0, numberTimesOptimal_ - 1), numberTimesOptimal_);
3322#endif
3323            currentDualTolerance_ = dualTolerance_ * pow(2.0, numberTimesOptimal_ - 1);
3324          }
3325          if (numberChanged_ > 0) {
3326            handler_->message(CLP_DUAL_ORIGINAL, messages_)
3327              << CoinMessageEol;
3328            //realDualInfeasibilities = sumDualInfeasibilities_;
3329          } else if (numberChangedBounds < -1) {
3330            // bounds were flipped
3331            //gutsOfSolution(2);
3332          }
3333          if (stateDualColumn_ >= 0)
3334            stateDualColumn_++;
3335#ifndef TRY_ABC_GUS
3336          int returnCode = bounceTolerances(1);
3337          if (returnCode)
3338#endif // always go to primal
3339            problemStatus_ = 10;
3340          lastCleaned_ = numberIterations_;
3341#if CLP_CAN_HAVE_ZERO_OBJ
3342        } else {
3343          // no cost - skip checks
3344          problemStatus_ = 0;
3345        }
3346#endif
3347      } else {
3348        problemStatus_ = 0; // optimal
3349        if (lastCleaned_ < numberIterations_ && numberChanged_) {
3350          handler_->message(CLP_SIMPLEX_GIVINGUP, messages_)
3351            << CoinMessageEol;
3352        }
3353      }
3354    } else if (numberChangedBounds < -1) {
3355      // some flipped
3356      bounceTolerances(2);
3357#if ABC_NORMAL_DEBUG > 0
3358      printf("numberChangedBounds %d numberAtFakeBound %d at line %d in file %s\n",
3359        numberChangedBounds, numberAtFakeBound(), __LINE__, __FILE__);
3360#endif
3361    } else if (numberChangedBounds > 0 || numberAtFakeBound()) {
3362      handler_->message(CLP_DUAL_BOUNDS, messages_)
3363        << currentDualBound_
3364        << CoinMessageEol;
3365#if ABC_NORMAL_DEBUG > 0
3366      printf("changing bounds\n");
3367      printf("numberChangedBounds %d numberAtFakeBound %d at line %d in file %s\n",
3368        numberChangedBounds, numberAtFakeBound(), __LINE__, __FILE__);
3369#endif
3370      if (!numberDualInfeasibilities_ && !numberPrimalInfeasibilities_) {
3371        // check unbounded
3372        // find a variable with bad dj
3373        int iChosen = -1;
3374        if ((specialOptions_ & 0x03000000) == 0) {
3375          double largest = 100.0 * primalTolerance_;
3376          for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
3377            double djValue = abcDj_[iSequence];
3378            double originalLo = lowerSaved_[iSequence];
3379            double originalUp = upperSaved_[iSequence];
3380            if (fabs(djValue) > fabs(largest)) {
3381              if (getInternalStatus(iSequence) != basic) {
3382                if (djValue > 0 && originalLo < -1.0e20) {
3383                  if (djValue > fabs(largest)) {
3384                    largest = djValue;
3385                    iChosen = iSequence;
3386                  }
3387                } else if (djValue < 0 && originalUp > 1.0e20) {
3388                  if (-djValue > fabs(largest)) {
3389                    largest = djValue;
3390                    iChosen = iSequence;
3391                  }
3392                }
3393              }
3394            }
3395          }
3396        }
3397        if (iChosen >= 0) {
3398          int iVector = getAvailableArray();
3399          unpack(usefulArray_[iVector], iChosen);
3400          //double changeCost;
3401          problemStatus_ = checkUnbounded(usefulArray_[iVector], 0.0 /*changeCost*/);
3402          usefulArray_[iVector].clear();
3403          setAvailableArray(iVector);
3404        }
3405      }
3406      gutsOfSolution(3);
3407      //bounceTolerances(2);
3408      //realDualInfeasibilities = sumDualInfeasibilities_;
3409    }
3410    // unflag all variables (we may want to wait a bit?)
3411    if (unflagVariables) {
3412      int numberFlagged = numberFlagged_;
3413      numberFlagged_ = 0;
3414      for (int iRow = 0; iRow < numberRows_; iRow++) {
3415        int iPivot = abcPivotVariable_[iRow];
3416        if (flagged(iPivot)) {
3417          clearFlagged(iPivot);
3418        }
3419      }
3420#if ABC_NORMAL_DEBUG > 3
3421      if (numberFlagged) {
3422        printf("unflagging %d variables - tentativeStatus %d probStat %d ninf %d nopt %d\n", numberFlagged, tentativeStatus,
3423          problemStatus_, numberPrimalInfeasibilities_,
3424          numberTimesOptimal_);
3425      }
3426#endif
3427      unflagVariables = numberFlagged > 0;
3428      if (numberFlagged && !numberPivots) {
3429        /* looks like trouble as we have not done any iterations.
3430           Try changing pivot tolerance then give it a few goes and give up */
3431        if (abcFactorization_->pivotTolerance() < 0.9) {
3432          abcFactorization_->pivotTolerance(0.99);
3433          problemStatus_ = -1;
3434        } else if (numberTimesOptimal_ < 3) {
3435          numberTimesOptimal_++;
3436          problemStatus_ = -1;
3437        } else {
3438          unflagVariables = false;
3439          //secondaryStatus_ = 1; // and say probably infeasible
3440          if ((moreSpecialOptions_ & 256) == 0 || (abcState_ & CLP_ABC_BEEN_FEASIBLE) != 0) {
3441            // try primal
3442            problemStatus_ = 10;
3443          } else {
3444            // almost certainly infeasible
3445            problemStatus_ = 1;
3446          }
3447#if ABC_NORMAL_DEBUG > 3
3448          printf("returning at %d\n", __LINE__);
3449#endif
3450        }
3451      }
3452    }
3453  }
3454  if (problemStatus_ < 0) {
3455    saveGoodStatus();
3456    if (weightsSaved) {
3457      // restore weights (if saved) - also recompute infeasibility list
3458      abcDualRowPivot_->saveWeights(this, weightsSaved);
3459    } else {
3460      abcDualRowPivot_->recomputeInfeasibilities();
3461    }
3462  }
3463  // see if cutoff reached
3464  checkCutoff(false);
3465  // If we are in trouble and in branch and bound give up
3466  if ((specialOptions_ & 1024) != 0) {
3467    int looksBad = 0;
3468    if (largestPrimalError_ * largestDualError_ > 1.0e2) {
3469      looksBad = 1;
3470    } else if (largestPrimalError_ > 1.0e-2
3471      && rawObjectiveValue_ > CoinMin(1.0e15, 1.0e3 * limit)) {
3472      looksBad = 2;
3473    }
3474    if (looksBad) {
3475      if (abcFactorization_->pivotTolerance() < 0.9) {
3476        // up tolerance
3477        abcFactorization_->pivotTolerance(CoinMin(abcFactorization_->pivotTolerance() * 1.05 + 0.02, 0.91));
3478      } else if (numberIterations_ > 10000) {
3479#if ABC_NORMAL_DEBUG > 0
3480        if (handler_->logLevel() > 2)
3481          printf("bad dual - saying infeasible %d\n", looksBad);
3482#endif
3483        problemStatus_ = 1;
3484        secondaryStatus_ = 1; // and say was on cutoff
3485      } else if (largestPrimalError_ > 1.0e5) {
3486#if ABC_NORMAL_DEBUG > 0
3487        if (handler_->logLevel() > 2)
3488          printf("bad dual - going to primal %d %g\n", looksBad, largestPrimalError_);
3489#endif
3490        allSlackBasis();
3491        problemStatus_ = 10;
3492      }
3493    }
3494  }
3495  //if (problemStatus_ < 0 && !changeMade_) {
3496  //problemStatus_ = 4; // unknown
3497  //}
3498  lastGoodIteration_ = numberIterations_;
3499  if (numberIterations_ > lastBadIteration_ + 200) {
3500    moreSpecialOptions_ &= ~16; // clear check accuracy flag
3501    if (numberFlagged_) {
3502      numberFlagged_ = 0;
3503      for (int iRow = 0; iRow < numberRows_; iRow++) {
3504        int iPivot = abcPivotVariable_[iRow];
3505        if (flagged(iPivot)) {
3506          clearFlagged(iPivot);
3507        }
3508      }
3509    }
3510  }
3511  if (problemStatus_ < 0 && ignoreStuff != 2) {
3512    //sumDualInfeasibilities_ = realDualInfeasibilities; // back to say be careful
3513    if (sumDualInfeasibilities_) {
3514      numberDualInfeasibilities_ = 1;
3515    } else if (!numberPrimalInfeasibilities_) {
3516      // Looks good
3517      problemStatus_ = 0;
3518    }
3519  }
3520  double thisObj = abcProgress_.lastObjective(0);
3521  double lastObj = abcProgress_.lastObjective(1);
3522  if (lastObj > thisObj + 1.0e-4 * CoinMax(fabs(thisObj), fabs(lastObj)) + 1.0e-4 && lastFirstFree_ < 0) {
3523    //printf("BAD - objective backwards\n");
3524    //bounceTolerances(3);
3525    numberBlackMarks += 3;
3526  }
3527  if (numberBlackMarks > 0 && numberIterations_ > baseIteration_) {
3528    // be very careful
3529    int maxFactor = abcFactorization_->maximumPivots();
3530    if (maxFactor > 10) {
3531      if ((stateOfProblem_ & PESSIMISTIC) == 0 || true) {
3532        if (largestDualError_ > 10.0 * CoinMax(lastDualError_, 1.0e-6)
3533          || largestPrimalError_ > 10.0 * CoinMax(lastPrimalError_, 1.0e-6))
3534          maxFactor = CoinMin(maxFactor, 20);
3535        forceFactorization_ = CoinMax(1, (maxFactor >> 1));
3536      } else if (numberBlackMarks > 2) {
3537        forceFactorization_ = CoinMax(1, (forceFactorization_ >> 1));
3538      }
3539    }
3540    stateOfProblem_ |= PESSIMISTIC;
3541    if (numberBlackMarks >= 10)
3542      forceFactorization_ = 1;
3543  }
3544  lastPrimalError_ = largestPrimalError_;
3545  lastDualError_ = largestDualError_;
3546  lastFirstFree_ = firstFree_;
3547  if (ignoreStuff == 2) {
3548    // in values pass
3549    lastFirstFree_ = -1;
3550    // firstFree_=-1;
3551  }
3552  if (alphaAccuracy_ > 0.0)
3553    alphaAccuracy_ = 1.0;
3554  // If we are stopping - use plausible objective
3555  // Maybe only in fast dual
3556  if (problemStatus_ > 2)
3557    rawObjectiveValue_ = approximateObjective;
3558  if (problemStatus_ == 1 && (progressFlag_ & 8) != 0 && fabs(rawObjectiveValue_) > 1.0e10)
3559    problemStatus_ = 10; // infeasible - but has looked feasible
3560  checkMoveBack(true);
3561#if 0
3562  {
3563    double * array = usefulArray_[arrayForTableauRow_].denseVector();
3564    for (int iPivot = 0; iPivot < numberRows_; iPivot++) {
3565      int iSequence = abcPivotVariable_[iPivot];
3566      unpack(usefulArray_[arrayForTableauRow_],iSequence);
3567      abcFactorization_->updateColumn(usefulArray_[arrayForTableauRow_]);
3568      assert (fabs(array[iPivot] - 1.0) < 1.0e-4);
3569      array[iPivot] = 0.0;
3570      for (int i = 0; i < numberRows_; i++)
3571        assert (fabs(array[i]) < 1.0e-4);
3572      usefulArray_[arrayForTableauRow_].clear();
3573    }
3574  }
3575#endif
3576#if 0
3577  {
3578    // save status, cost and solution
3579    unsigned char * saveStatus = CoinCopyOfArray(internalStatus_,numberTotal_);
3580    double * saveSolution = CoinCopyOfArray(abcSolution_,numberTotal_);
3581    double * saveCost = CoinCopyOfArray(abcCost_,numberTotal_);
3582    printf("TestA sumDual %g sumPrimal %g\n",sumDualInfeasibilities_,sumPrimalInfeasibilities_);
3583    if (perturbation_>100)
3584      CoinAbcMemcpy(abcCost_,costSaved_,numberTotal_);
3585    else
3586      CoinAbcMemcpy(abcCost_,abcPerturbation_,numberTotal_);
3587    moveToBasic();
3588    gutsOfSolution(3);
3589    for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
3590      bool changed=fabs(saveCost[iSequence]-abcCost_[iSequence])>1.0e-10;
3591      bool bad=false;
3592      int iStatus=internalStatus_[iSequence]&7;
3593      if (getInternalStatus(iSequence) != basic && !flagged(iSequence)) {
3594        // not basic
3595        double distanceUp = abcUpper_[iSequence] - abcSolution_[iSequence];
3596        double distanceDown = abcSolution_[iSequence] - abcLower_[iSequence];
3597        double value = abcDj_[iSequence];
3598        if (distanceUp > primalTolerance_) {
3599          // Check if "free"
3600          if (distanceDown <= primalTolerance_) {
3601            // should not be negative
3602            if (value < 0.0) {
3603              value = - value;
3604              if (value > currentDualTolerance_) {
3605                bad=true;
3606              }
3607            }
3608          } else {
3609            // free
3610            value=fabs(value);
3611            if (value > currentDualTolerance_) {
3612              bad=true;
3613            }
3614          }
3615        } else if (distanceDown > primalTolerance_) {
3616          // should not be positive
3617          if (value > 0.0) {
3618            if (value > currentDualTolerance_) {
3619              bad=true;
3620            }
3621          }
3622        }
3623      }
3624      if (changed||bad)
3625        printf("seq %d status %d current cost %g original %g dj %g %s\n",
3626               iSequence,iStatus,saveCost[iSequence],abcCost_[iSequence],
3627               abcDj_[iSequence],bad ? "(bad)" : "");
3628    }
3629    printf("TestB sumDual %g sumPrimal %g\n",sumDualInfeasibilities_,sumPrimalInfeasibilities_);
3630    if (numberDualInfeasibilities_) {
3631      makeNonFreeVariablesDualFeasible();
3632      moveToBasic();
3633      gutsOfSolution(3);
3634      printf("TestC sumDual %g sumPrimal %g\n",sumDualInfeasibilities_,sumPrimalInfeasibilities_);
3635    }
3636    CoinAbcMemcpy(internalStatus_,saveStatus,numberTotal_);
3637    CoinAbcMemcpy(abcSolution_,saveSolution,numberTotal_);
3638    CoinAbcMemcpy(abcCost_,saveCost,numberTotal_);
3639    delete [] saveStatus;
3640    delete [] saveSolution;
3641    delete [] saveCost;
3642    moveToBasic();
3643    gutsOfSolution(3);
3644    printf("TestD sumDual %g sumPrimal %g\n",sumDualInfeasibilities_,sumPrimalInfeasibilities_);
3645  }
3646#endif
3647#if 0
3648  {
3649    const double multiplier[] = { 1.0, -1.0}; // sign?????????????
3650    double dualT = - currentDualTolerance_;
3651    double largestCostDifference=0.1*currentDualTolerance_;
3652    int numberCostDifference=0;
3653    double largestBadDj=0.0;
3654    int numberBadDj=0;
3655    double sumBadDj=0.0;
3656    double sumCanMoveBad=0.0;
3657    int numberCanMoveBad=0;
3658    double sumCanMoveOthers=0.0;
3659    int numberCanMoveOthers=0;
3660    for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
3661      double costDifference=abcPerturbation_[iSequence]-abcCost_[iSequence];
3662      if (fabs(costDifference)>0.1*currentDualTolerance_) {
3663        if (fabs(costDifference)>largestCostDifference)
3664          largestCostDifference=fabs(costDifference);
3665        numberCostDifference++;
3666      }
3667      unsigned char iStatus = internalStatus_[iSequence] & 7;
3668      if (iStatus<4) {
3669        double mult = multiplier[iStatus];
3670        double dj = abcDj_[iSequence] * mult;
3671        if (dj < dualT) {
3672          sumBadDj+=dualT-dj;
3673          numberBadDj++;
3674          if (dualT-dj>largestBadDj)
3675            largestBadDj=dualT-dj;
3676          if (costDifference*mult>0.0) {
3677            sumCanMoveBad+=costDifference*mult;
3678            numberCanMoveBad++;
3679            //dj+=costDifference*mult;
3680            abcDj_[iSequence]+=costDifference;
3681            abcCost_[iSequence]=abcPerturbation_[iSequence];
3682          }
3683        } else if (costDifference*mult>0.0) {
3684          sumCanMoveOthers+=costDifference*mult;
3685          numberCanMoveOthers++;
3686          //dj+=costDifference*mult;
3687          abcDj_[iSequence]+=costDifference;
3688          abcCost_[iSequence]=abcPerturbation_[iSequence];
3689        }
3690      }
3691    }
3692    if (numberCostDifference+numberBadDj+numberCanMoveBad+numberCanMoveOthers) {
3693      printf("largest diff %g (%d), largest bad dj %g (%d - sum %g), can move bad %g (%d), good %g (%d) - %d at fake bound\n",largestCostDifference,numberCostDifference,largestBadDj,
3694             numberBadDj, sumBadDj,sumCanMoveBad=0.0,numberCanMoveBad,
3695             sumCanMoveOthers,numberCanMoveOthers,numberAtFakeBound());
3696    }
3697  }
3698#endif
3699}
3700// see if cutoff reached
3701bool AbcSimplexDual::checkCutoff(bool computeObjective)
3702{
3703  double objValue = COIN_DBL_MAX;
3704  if (computeObjective)
3705    objValue = computeInternalObjectiveValue();
3706  else
3707    objValue = minimizationObjectiveValue();
3708  double limit = 0.0;
3709  getDblParam(ClpDualObjectiveLimit, limit);
3710  // Looks as if numberAtFakeBound() being computed too often
3711  if (fabs(limit) < 1.0e30 && objValue > limit && !numberAtFakeBound()) {
3712    bool looksInfeasible = !numberDualInfeasibilities_;
3713    if (objValue > limit + fabs(0.1 * limit) + 1.0e2 * sumDualInfeasibilities_ + 1.0e4 && sumDualInfeasibilities_ < largestDualError_ && numberIterations_ > 0.5 * numberRows_ + 1000)
3714      looksInfeasible = true;
3715    if (looksInfeasible && !computeObjective)
3716      objValue = computeInternalObjectiveValue();
3717    if (looksInfeasible) {
3718      // Even if not perturbed internal costs may have changed
3719      // be careful
3720      if (objValue > limit) {
3721        problemStatus_ = 1;
3722        secondaryStatus_ = 1; // and say was on cutoff
3723        return true;
3724      }
3725    }
3726  }
3727  return false;
3728}
3729#define CLEANUP_DJS 0
3730/* While updateDualsInDual sees what effect is of flip
3731   this does actual flipping.
3732   If change >0.0 then value in array >0.0 => from lower to upper
3733   returns 3 if skip this iteration and re-factorize
3734*/
3735int AbcSimplexDual::flipBounds()
3736{
3737  CoinIndexedVector &array = usefulArray_[arrayForFlipBounds_];
3738  CoinIndexedVector &output = usefulArray_[arrayForFlipRhs_];
3739  output.clear();
3740  int numberFlipped = array.getNumElements();
3741  if (numberFlipped) {
3742#if CLP_CAN_HAVE_ZERO_OBJ > 1
3743    if ((specialOptions_ & 2097152) == 0) {
3744#endif
3745      // do actual flips
3746      // do in parallel once got spare space in factorization
3747#if ABC_PARALLEL == 2
3748#endif
3749      int number = array.getNumElements();
3750      const int *COIN_RESTRICT which = array.getIndices();
3751      double *COIN_RESTRICT work = array.denseVector();
3752
3753      double *COIN_RESTRICT solution = abcSolution_;
3754      const double *COIN_RESTRICT lower = abcLower_;
3755      const double *COIN_RESTRICT upper = abcUpper_;
3756#if CLEANUP_DJS
3757      double *COIN_RESTRICT cost = abcCost_;
3758      double *COIN_RESTRICT dj = abcDj_;
3759      // see if we can clean up djs
3760      // may have perturbation
3761      const double *COIN_RESTRICT originalCost = abcPerturbation_;
3762#else
3763    const double *COIN_RESTRICT cost = abcCost_;
3764    const double *COIN_RESTRICT dj = abcDj_;
3765#endif
3766      const double multiplier[] = { 1.0, -1.0 };
3767      for (int i = 0; i < number; i++) {
3768        double value = work[i] * theta_;
3769        work[i] = 0.0;
3770        int iSequence = which[i];
3771        unsigned char iStatus = internalStatus_[iSequence] & 3;
3772        assert((internalStatus_[iSequence] & 7) == 0 || (internalStatus_[iSequence] & 7) == 1);
3773        double mult = multiplier[iStatus];
3774        double djValue = dj[iSequence] * mult - value;
3775        //double perturbedTolerance = abcPerturbation_[iSequence];
3776        if (djValue < -currentDualTolerance_) {
3777          // might just happen - if so just skip
3778          assert(iSequence != sequenceIn_);
3779          double movement = (upper[iSequence] - lower[iSequence]) * mult;
3780          if (iStatus == 0) {
3781            // at lower bound
3782            // to upper bound
3783            setInternalStatus(iSequence, atUpperBound);
3784            solution[iSequence] = upper[iSequence];
3785#if CLEANUP_DJS
3786            if (cost[iSequence] < originalCost[iSequence]) {
3787              double difference = CoinMax(cost[iSequence] - originalCost[iSequence], djValue);
3788              dj[iSequence] -= difference;
3789              cost[iSequence] -= difference;
3790            }
3791#endif
3792          } else {
3793            // at upper bound
3794            // to lower bound
3795            setInternalStatus(iSequence, atLowerBound);
3796            solution[iSequence] = lower[iSequence];
3797#if CLEANUP_DJS
3798            if (cost[iSequence] > originalCost[iSequence]) {
3799              double difference = CoinMin(cost[iSequence] - originalCost[iSequence], -djValue);
3800              dj[iSequence] -= difference;
3801              cost[iSequence] -= difference;
3802            }
3803#endif
3804          }
3805          abcMatrix_->add(output, iSequence, movement);
3806          objectiveChange_ += cost[iSequence] * movement;
3807        }
3808      }
3809      array.setNumElements(0);
3810      // double check something flipped
3811      if (output.getNumElements()) {
3812#if ABC_PARALLEL == 0
3813        abcFactorization_->updateColumn(output);
3814#else
3815      assert(FACTOR_CPU > 2);
3816      abcFactorization_->updateColumnCpu(output, 1);
3817#endif
3818        abcDualRowPivot_->updatePrimalSolution(output, 1.0);
3819      }
3820#if CLP_CAN_HAVE_ZERO_OBJ > 1
3821    } else {
3822      clearArrays(3);
3823    }
3824#endif
3825  }
3826  // recompute dualOut_
3827  valueOut_ = solutionBasic_[pivotRow_];
3828  if (directionOut_ < 0) {
3829    dualOut_ = valueOut_ - upperOut_;
3830  } else {
3831    dualOut_ = lowerOut_ - valueOut_;
3832  }
3833#if 0
3834  // amount primal will move
3835  movement_ = -dualOut_ * directionOut_ / alpha_;
3836  double movementOld = oldDualOut * directionOut_ / alpha_;
3837  // so objective should increase by fabs(dj)*movement_
3838  // but we already have objective change - so check will be good
3839  if (objectiveChange_ + fabs(movementOld * dualIn_) < -CoinMax(1.0e-5, 1.0e-12 * fabs(rawObjectiveValue_))) {
3840#if ABC_NORMAL_DEBUG > 3
3841    if (handler_->logLevel() & 32)
3842      printf("movement %g, movementOld %g swap change %g, rest %g  * %g\n",
3843             objectiveChange_ + fabs(movement_ * dualIn_),
3844             movementOld,objectiveChange_, movement_, dualIn_);
3845#endif
3846  }
3847  if (0) {
3848    if(abcFactorization_->pivots()) {
3849      // going backwards - factorize
3850      problemStatus_ = -2; // factorize now
3851      stateOfIteration_=2;
3852    }
3853  }
3854#endif
3855  return numberFlipped;
3856}
3857/* Undo a flip
3858 */
3859void AbcSimplexDual::flipBack(int number)
3860{
3861  CoinIndexedVector &array = usefulArray_[arrayForFlipBounds_];
3862  const int *COIN_RESTRICT which = array.getIndices();
3863  double *COIN_RESTRICT solution = abcSolution_;
3864  const double *COIN_RESTRICT lower = abcLower_;
3865  const double *COIN_RESTRICT upper = abcUpper_;
3866  for (int i = 0; i < number; i++) {
3867    int iSequence = which[i];
3868    unsigned char iStatus = internalStatus_[iSequence] & 3;
3869    assert((internalStatus_[iSequence] & 7) == 0 || (internalStatus_[iSequence] & 7) == 1);
3870    if (iStatus == 0) {
3871      // at lower bound
3872      // to upper bound
3873      setInternalStatus(iSequence, atUpperBound);
3874      solution[iSequence] = upper[iSequence];
3875    } else {
3876      // at upper bound
3877      // to lower bound
3878      setInternalStatus(iSequence, atLowerBound);
3879      solution[iSequence] = lower[iSequence];
3880    }
3881  }
3882}
3883// Restores bound to original bound
3884void AbcSimplexDual::originalBound(int iSequence)
3885{
3886  if (getFakeBound(iSequence) != noFake) {
3887    numberFake_--;
3888    ;
3889    setFakeBound(iSequence, noFake);
3890    abcLower_[iSequence] = lowerSaved_[iSequence];
3891    abcUpper_[iSequence] = upperSaved_[iSequence];
3892  }
3893}
3894/* As changeBounds but just changes new bounds for a single variable.
3895   Returns true if change */
3896bool AbcSimplexDual::changeBound(int iSequence)
3897{
3898  // old values
3899  double oldLower = abcLower_[iSequence];
3900  double oldUpper = abcUpper_[iSequence];
3901  double value = abcSolution_[iSequence];
3902  bool modified = false;
3903  // original values
3904  double lowerValue = lowerSaved_[iSequence];
3905  double upperValue = upperSaved_[iSequence];
3906  setFakeBound(iSequence, noFake);
3907  if (value == oldLower) {
3908    if (upperValue > oldLower + currentDualBound_) {
3909      abcUpper_[iSequence] = oldLower + currentDualBound_;
3910      setFakeBound(iSequence, upperFake);
3911      assert(abcLower_[iSequence] == lowerSaved_[iSequence]);
3912      assert(abcUpper_[iSequence] < upperSaved_[iSequence]);
3913      modified = true;
3914      numberFake_++;
3915    }
3916  } else if (value == oldUpper) {
3917    if (lowerValue < oldUpper - currentDualBound_) {
3918      abcLower_[iSequence] = oldUpper - currentDualBound_;
3919      setFakeBound(iSequence, lowerFake);
3920      assert(abcLower_[iSequence] > lowerSaved_[iSequence]);
3921      assert(abcUpper_[iSequence] == upperSaved_[iSequence]);
3922      modified = true;
3923      numberFake_++;
3924    }
3925  } else {
3926    assert(value == oldLower || value == oldUpper);
3927  }
3928  return modified;
3929}
3930/* Checks number of variables at fake bounds.  This is used by fastDual
3931   so can exit gracefully before end */
3932int AbcSimplexDual::numberAtFakeBound()
3933{
3934  int numberFake = 0;
3935
3936  if (numberFake_) {
3937    for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
3938      FakeBound bound = getFakeBound(iSequence);
3939      switch (getInternalStatus(iSequence)) {
3940
3941      case basic:
3942        break;
3943      case isFree:
3944      case superBasic:
3945      case AbcSimplex::isFixed:
3946        //setFakeBound (iSequence, noFake);
3947        break;
3948      case atUpperBound:
3949        if (bound == upperFake || bound == bothFake)
3950          numberFake++;
3951        break;
3952      case atLowerBound:
3953        if (bound == lowerFake || bound == bothFake)
3954          numberFake++;
3955        break;
3956      }
3957    }
3958  }
3959  return numberFake;
3960}
3961int AbcSimplexDual::resetFakeBounds(int type)
3962{
3963  if (type == 0 || type == 1) {
3964    // put back original bounds and then check
3965    //CoinAbcMemcpy(abcLower_,lowerSaved_,numberTotal_);
3966    //CoinAbcMemcpy(abcUpper_,upperSaved_,numberTotal_);
3967    double dummyChangeCost = 0.0;
3968    return changeBounds(3 + type, dummyChangeCost);
3969  } else if (type < 0) {
3970#ifndef NDEBUG
3971    // just check
3972    int iSequence;
3973    int nFake = 0;
3974    int nErrors = 0;
3975    int nSuperBasic = 0;
3976    int nWarnings = 0;
3977    for (iSequence = 0; iSequence < numberTotal_; iSequence++) {
3978      FakeBound fakeStatus = getFakeBound(iSequence);
3979      Status status = getInternalStatus(iSequence);
3980      bool isFake = false;
3981#ifdef CLP_INVESTIGATE
3982      char RC = 'C';
3983      int jSequence = iSequence;
3984      if (jSequence >= numberColumns_) {
3985        RC = 'R';
3986        jSequence -= numberColumns_;
3987      }
3988#endif
3989      double lowerValue = lowerSaved_[iSequence];
3990      double upperValue = upperSaved_[iSequence];
3991      double value = abcSolution_[iSequence];
3992      CoinRelFltEq equal;
3993      if (status == atUpperBound || status == atLowerBound) {
3994        if (fakeStatus == AbcSimplexDual::upperFake) {
3995          if (!equal(abcUpper_[iSequence], (lowerValue + currentDualBound_)) || !(equal(abcUpper_[iSequence], value) || equal(abcLower_[iSequence], value))) {
3996            nErrors++;
3997#ifdef CLP_INVESTIGATE
3998            printf("** upperFake %c%d %g <= %g <= %g true %g, %g\n",
3999              RC, jSequence, abcLower_[iSequence], abcSolution_[iSequence],
4000              abcUpper_[iSequence], lowerValue, upperValue);
4001#endif
4002          }
4003          isFake = true;
4004          ;
4005        } else if (fakeStatus == AbcSimplexDual::lowerFake) {
4006          if (!equal(abcLower_[iSequence], (upperValue - currentDualBound_)) || !(equal(abcUpper_[iSequence], value) || equal(abcLower_[iSequence], value))) {
4007            nErrors++;
4008#ifdef CLP_INVESTIGATE
4009            printf("** lowerFake %c%d %g <= %g <= %g true %g, %g\n",
4010              RC, jSequence, abcLower_[iSequence], abcSolution_[iSequence],
4011              abcUpper_[iSequence], lowerValue, upperValue);
4012#endif
4013          }
4014          isFake = true;
4015          ;
4016        } else if (fakeStatus == AbcSimplexDual::bothFake) {
4017          nWarnings++;
4018#ifdef CLP_INVESTIGATE
4019          printf("** %d at bothFake?\n", iSequence);
4020#endif
4021        } else if (abcUpper_[iSequence] - abcLower_[iSequence] > 2.0 * currentDualBound_) {
4022          nErrors++;
4023#ifdef CLP_INVESTIGATE
4024          printf("** noFake! %c%d %g <= %g <= %g true %g, %g\n",
4025            RC, jSequence, abcLower_[iSequence], abcSolution_[iSequence],
4026            abcUpper_[iSequence], lowerValue, upperValue);
4027#endif
4028        }
4029      } else if (status == superBasic || status == isFree) {
4030        nSuperBasic++;
4031        //printf("** free or superbasic %c%d %g <= %g <= %g true %g, %g - status %d\n",
4032        //     RC,jSequence,abcLower_[iSequence],abcSolution_[iSequence],
4033        //     abcUpper_[iSequence],lowerValue,upperValue,status);
4034      } else if (status == basic) {
4035        bool odd = false;
4036        if (!equal(abcLower_[iSequence], lowerValue))
4037          odd = true;
4038        if (!equal(abcUpper_[iSequence], upperValue))
4039          odd = true;
4040        if (odd) {
4041#ifdef CLP_INVESTIGATE
4042          printf("** basic %c%d %g <= %g <= %g true %g, %g\n",
4043            RC, jSequence, abcLower_[iSequence], abcSolution_[iSequence],
4044            abcUpper_[iSequence], lowerValue, upperValue);
4045#endif
4046          nWarnings++;
4047        }
4048      } else if (status == isFixed) {
4049        if (!equal(abcUpper_[iSequence], abcLower_[iSequence])) {
4050          nErrors++;
4051#ifdef CLP_INVESTIGATE
4052          printf("** fixed! %c%d %g <= %g <= %g true %g, %g\n",
4053            RC, jSequence, abcLower_[iSequence], abcSolution_[iSequence],
4054            abcUpper_[iSequence], lowerValue, upperValue);
4055#endif
4056        }
4057      }
4058      if (isFake) {
4059        nFake++;
4060      } else {
4061        if (fakeStatus != AbcSimplexDual::noFake) {
4062          nErrors++;
4063#ifdef CLP_INVESTIGATE
4064          printf("** bad fake status %c%d %d\n",
4065            RC, jSequence, fakeStatus);
4066#endif
4067        }
4068      }
4069    }
4070    if (nFake != numberFake_) {
4071#ifdef CLP_INVESTIGATE
4072      printf("nfake %d numberFake %d\n", nFake, numberFake_);
4073#endif
4074      nErrors++;
4075    }
4076    if (nErrors || type <= -1000) {
4077#ifdef CLP_INVESTIGATE
4078      printf("%d errors, %d warnings, %d free/superbasic, %d fake\n",
4079        nErrors, nWarnings, nSuperBasic, numberFake_);
4080      printf("currentDualBound %g\n",
4081        currentDualBound_);
4082#endif
4083      if (type <= -1000) {
4084        iSequence = -type;
4085        iSequence -= 1000;
4086#ifdef CLP_INVESTIGATE
4087        char RC = 'C';
4088        int jSequence = iSequence;
4089        if (jSequence >= numberColumns_) {
4090          RC = 'R';
4091          jSequence -= numberColumns_;
4092        }
4093        double lowerValue = lowerSaved_[iSequence];
4094        double upperValue = upperSaved_[iSequence];
4095        printf("*** movement>1.0e30 for  %c%d %g <= %g <= %g true %g, %g - status %d\n",
4096          RC, jSequence, abcLower_[iSequence], abcSolution_[iSequence],
4097          abcUpper_[iSequence], lowerValue, upperValue, internalStatus_[iSequence]);
4098#endif
4099        assert(nErrors); // should have been picked up
4100      }
4101      assert(!nErrors);
4102    }
4103#endif
4104  } else {
4105    CoinAbcMemcpy(abcLower_, lowerSaved_, numberTotal_);
4106    CoinAbcMemcpy(abcUpper_, upperSaved_, numberTotal_);
4107    // reset using status
4108    numberFake_ = 0;
4109    for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
4110      FakeBound fakeStatus = getFakeBound(iSequence);
4111      if (fakeStatus != AbcSimplexDual::noFake) {
4112        Status status = getInternalStatus(iSequence);
4113        if (status == basic) {
4114          setFakeBound(iSequence, AbcSimplexDual::noFake);
4115          continue;
4116        }
4117        double lowerValue = abcLower_[iSequence];
4118        double upperValue = abcUpper_[iSequence];
4119        double value = abcSolution_[iSequence];
4120        numberFake_++;
4121        if (fakeStatus == AbcSimplexDual::upperFake) {
4122          abcUpper_[iSequence] = lowerValue + currentDualBound_;
4123          if (status == AbcSimplex::atLowerBound) {
4124            abcSolution_[iSequence] = lowerValue;
4125          } else if (status == AbcSimplex::atUpperBound) {
4126            abcSolution_[iSequence] = abcUpper_[iSequence];
4127          } else {
4128            abort();
4129          }
4130        } else if (fakeStatus == AbcSimplexDual::lowerFake) {
4131          abcLower_[iSequence] = upperValue - currentDualBound_;
4132          if (status == AbcSimplex::atLowerBound) {
4133            abcSolution_[iSequence] = abcLower_[iSequence];
4134          } else if (status == AbcSimplex::atUpperBound) {
4135            abcSolution_[iSequence] = upperValue;
4136          } else {
4137            abort();
4138          }
4139        } else {
4140          assert(fakeStatus == AbcSimplexDual::bothFake);
4141          if (status == AbcSimplex::atLowerBound) {
4142            abcLower_[iSequence] = value;
4143            abcUpper_[iSequence] = value + currentDualBound_;
4144          } else if (status == AbcSimplex::atUpperBound) {
4145            abcUpper_[iSequence] = value;
4146            abcLower_[iSequence] = value - currentDualBound_;
4147          } else if (status == AbcSimplex::isFree || status == AbcSimplex::superBasic) {
4148            abcLower_[iSequence] = value - 0.5 * currentDualBound_;
4149            abcUpper_[iSequence] = value + 0.5 * currentDualBound_;
4150          } else {
4151            abort();
4152          }
4153        }
4154      }
4155    }
4156  }
4157  return -1;
4158}
4159// maybe could JUST check against perturbationArray
4160void AbcSimplex::checkDualSolutionPlusFake()
4161{
4162
4163  sumDualInfeasibilities_ = 0.0;
4164  numberDualInfeasibilities_ = 0;
4165  sumFakeInfeasibilities_ = 0.0;
4166  int firstFreePrimal = -1;
4167  int firstFreeDual = -1;
4168  int numberSuperBasicWithDj = 0;
4169  bestPossibleImprovement_ = 0.0;
4170  // we can't really trust infeasibilities if there is dual error
4171  double error = CoinMin(1.0e-2, largestDualError_);
4172  // allow tolerance at least slightly bigger than standard
4173  double relaxedTolerance = currentDualTolerance_ + error;
4174  // allow bigger tolerance for possible improvement
4175  double possTolerance = 5.0 * relaxedTolerance;
4176  sumOfRelaxedDualInfeasibilities_ = 0.0;
4177  //rawObjectiveValue_ -=sumNonBasicCosts_;
4178  sumNonBasicCosts_ = 0.0;
4179  const double *COIN_RESTRICT abcLower = abcLower_;
4180  const double *COIN_RESTRICT abcUpper = abcUpper_;
4181  const double *COIN_RESTRICT abcSolution = abcSolution_;
4182#ifdef CLEANUP_DJS
4183  double *COIN_RESTRICT abcDj = abcDj_;
4184  double *COIN_RESTRICT abcCost = abcCost_;
4185  // see if we can clean up djs
4186  // may have perturbation
4187  const double *COIN_RESTRICT originalCost = abcPerturbation_;
4188#else
4189  const double *COIN_RESTRICT abcDj = abcDj_;
4190  const double *COIN_RESTRICT abcCost = abcCost_;
4191#endif
4192  //const double *  COIN_RESTRICT abcPerturbation = abcPerturbation_;
4193  int badDjSequence = -1;
4194  for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
4195    if (getInternalStatus(iSequence) != basic) {
4196      sumNonBasicCosts_ += abcSolution[iSequence] * abcCost[iSequence];
4197      // not basic
4198      double distanceUp = abcUpper[iSequence] - abcSolution[iSequence];
4199      double distanceDown = abcSolution[iSequence] - abcLower[iSequence];
4200      double value = abcDj[iSequence];
4201      if (distanceUp > primalTolerance_) {
4202        // Check if "free"
4203        if (distanceDown <= primalTolerance_) {
4204          // should not be negative
4205          if (value < 0.0) {
4206            value = -value;
4207            if (value > currentDualTolerance_) {
4208              bool treatAsSuperBasic = false;
4209              if (getInternalStatus(iSequence) == superBasic) {
4210#ifdef PAN
4211                if (fakeSuperBasic(iSequence) >= 0) {
4212#endif
4213                  treatAsSuperBasic = true;
4214                  numberSuperBasicWithDj++;
4215                  if (firstFreeDual < 0)
4216                    firstFreeDual = iSequence;
4217#ifdef PAN
4218                }
4219#endif
4220              }
4221              if (!treatAsSuperBasic) {
4222                sumDualInfeasibilities_ += value - currentDualTolerance_;
4223                badDjSequence = iSequence;
4224                if (value > possTolerance)
4225                  bestPossibleImprovement_ += CoinMin(distanceUp, 1.0e10) * value;
4226                if (value > relaxedTolerance)
4227                  sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
4228                numberDualInfeasibilities_++;
4229              }
4230            }
4231#ifdef CLEANUP_DJS
4232          } else {
4233            // adjust cost etc
4234            if (abcCost[iSequence] > originalCost[iSequence]) {
4235              double difference = CoinMin(abcCost[iSequence] - originalCost[iSequence], value);
4236              //printf("Lb %d changing cost from %g to %g, dj from %g to %g\n",
4237              //     iSequence,abcCost[iSequence],abcCost[iSequence]-difference,
4238              //     abcDj[iSequence],abcDj[iSequence]-difference);
4239              abcDj[iSequence] = value - difference;
4240              abcCost[iSequence] -= difference;
4241            }
4242#endif
4243          }
4244        } else {
4245          // free
4246          value = fabs(value);
4247          if (value > 1.0 * relaxedTolerance) {
4248            numberSuperBasicWithDj++;
4249            if (firstFreeDual < 0)
4250              firstFreeDual = iSequence;
4251          }
4252          if (value > currentDualTolerance_) {
4253            sumDualInfeasibilities_ += value - currentDualTolerance_;
4254            if (value > possTolerance)
4255              bestPossibleImprovement_ += CoinMin(distanceUp, 1.0e10) * value;
4256            if (value > relaxedTolerance)
4257              sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
4258            numberDualInfeasibilities_++;
4259          }
4260          if (firstFreePrimal < 0)
4261            firstFreePrimal = iSequence;
4262        }
4263      } else if (distanceDown > primalTolerance_) {
4264        // should not be positive
4265        if (value > 0.0) {
4266          if (value > currentDualTolerance_) {
4267            bool treatAsSuperBasic = false;
4268            if (getInternalStatus(iSequence) == superBasic) {
4269#ifdef PAN
4270              if (fakeSuperBasic(iSequence) >= 0) {
4271#endif
4272                treatAsSuperBasic = true;
4273                numberSuperBasicWithDj++;
4274                if (firstFreeDual < 0)
4275                  firstFreeDual = iSequence;
4276#ifdef PAN
4277              }
4278#endif
4279            }
4280            if (!treatAsSuperBasic) {
4281              sumDualInfeasibilities_ += value - currentDualTolerance_;
4282              badDjSequence = iSequence;
4283              if (value > possTolerance)
4284                bestPossibleImprovement_ += value * CoinMin(distanceDown, 1.0e10);
4285              if (value > relaxedTolerance)
4286                sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
4287              numberDualInfeasibilities_++;
4288            }
4289          }
4290#ifdef CLEANUP_DJS
4291        } else {
4292          // adjust cost etc
4293          if (abcCost[iSequence] < originalCost[iSequence]) {
4294            double difference = CoinMax(abcCost[iSequence] - originalCost[iSequence], value);
4295            //printf("Ub %d changing cost from %g to %g, dj from %g to %g\n",
4296            //       iSequence,abcCost[iSequence],abcCost[iSequence]-difference,
4297            //       abcDj[iSequence],abcDj[iSequence]-difference);
4298            abcDj[iSequence] = value - difference;
4299            abcCost[iSequence] -= difference;
4300          }
4301#endif
4302        }
4303      }
4304    }
4305  }
4306#ifdef CLEANUP_DJS
4307  if (algorithm_ < 0 && numberDualInfeasibilities_ == 1 && firstFreeDual < 0 && sumDualInfeasibilities_ < 1.0e-1) {
4308    // clean bad one
4309    assert(badDjSequence >= 0);
4310    abcCost[badDjSequence] -= abcDj_[badDjSequence];
4311    abcDj_[badDjSequence] = 0.0;
4312  }
4313#endif
4314  numberDualInfeasibilitiesWithoutFree_ = numberDualInfeasibilities_ - numberSuperBasicWithDj;
4315  numberFreeNonBasic_ = numberSuperBasicWithDj;
4316  firstFree_ = -1;
4317  if (algorithm_ < 0 && firstFreeDual >= 0) {
4318    // dual
4319    firstFree_ = firstFreeDual;
4320  } else if ((numberSuperBasicWithDj || (abcProgress_.lastIterationNumber(0) <= 0)) && algorithm_ > 0) {
4321    firstFree_ = firstFreePrimal;
4322  }
4323}
4324// Perturbs problem
4325void AbcSimplexDual::perturb(double factor)
4326{
4327  const double *COIN_RESTRICT lower = lowerSaved_;
4328  const double *COIN_RESTRICT upper = upperSaved_;
4329  const double *COIN_RESTRICT perturbation = abcPerturbation_;
4330  double *COIN_RESTRICT cost = abcCost_;
4331  for (int i = 0; i < numberTotal_; i++) {
4332    if (getInternalStatus(i) != basic) {
4333      if (fabs(lower[i]) < fabs(upper[i])) {
4334        cost[i] += factor * perturbation[i];
4335      } else if (upper[i] != COIN_DBL_MAX && upper[i] != lower[i]) {
4336        cost[i] -= factor * perturbation[i];
4337      } else {
4338        // free
4339        //cost[i]=costSaved[i];
4340      }
4341    }
4342  }
4343#ifdef CLEANUP_DJS
4344  // copy cost to perturbation
4345  CoinAbcMemcpy(abcPerturbation_, abcCost_, numberTotal_);
4346#endif
4347}
4348#if ABC_NORMAL_DEBUG > 0
4349#define PERT_STATISTICS
4350#endif
4351#ifdef PERT_STATISTICS
4352static void breakdown(const char *name, int numberLook, const double *region)
4353{
4354  double range[] = {
4355    -COIN_DBL_MAX,
4356    -1.0e15, -1.0e11, -1.0e8, -1.0e5, -3.0e4, -1.0e4, -3.0e3, -1.0e3, -3.0e2, -1.0e2, -3.0e1, -1.0e1, -3.0,
4357    -1.0, -3.0e-1,
4358    -1.0e-1, -3.0e-1, -1.0e-2, -3.0e-2, -1.0e-3, -3.0e-3, -1.0e-4, -3.0e-4, -1.0e-5, -1.0e-8, -1.0e-11, -1.0e-15,
4359    0.0,
4360    1.0e-15, 1.0e-11, 1.0e-8, 1.0e-5, 3.0e-5, 1.0e-4, 3.0e-4, 1.0e-3, 3.0e-3, 1.0e-2, 3.0e-2, 1.0e-1, 3.0e-1,
4361    1.0, 3.0,
4362    1.0e1, 3.0e1, 1.0e2, 3.0e2, 1.0e3, 3.0e3, 1.0e4, 3.0e4, 1.0e5, 1.0e8, 1.0e11, 1.0e15,
4363    COIN_DBL_MAX
4364  };
4365  int nRanges = static_cast< int >(sizeof(range) / sizeof(double));
4366  int *number = new int[nRanges];
4367  memset(number, 0, nRanges * sizeof(int));
4368  int *numberExact = new int[nRanges];
4369  memset(numberExact, 0, nRanges * sizeof(int));
4370  int i;
4371  for (i = 0; i < numberLook; i++) {
4372    double value = region[i];
4373    for (int j = 0; j < nRanges; j++) {
4374      if (value == range[j]) {
4375        numberExact[j]++;
4376        break;
4377      } else if (value < range[j]) {
4378        number[j]++;
4379        break;
4380      }
4381    }
4382  }
4383  printf("\n%s has %d entries\n", name, numberLook);
4384  for (i = 0; i < nRanges; i++) {
4385    if (number[i])
4386      printf("%d between %g and %g", number[i], range[i - 1], range[i]);
4387    if (numberExact[i]) {
4388      if (number[i])
4389        printf(", ");
4390      printf("%d exactly at %g", numberExact[i], range[i]);
4391    }
4392    if (number[i] + numberExact[i])
4393      printf("\n");
4394  }
4395  delete[] number;
4396  delete[] numberExact;
4397}
4398#endif
4399// Perturbs problem
4400void AbcSimplexDual::perturbB(double /*factorIn*/, int /*type*/)
4401{
4402  double overallMultiplier = (perturbation_ == 50 || perturbation_ > 54) ? 1.0 : 0.1;
4403  // dual perturbation
4404  double perturbation = 1.0e-20;
4405  // maximum fraction of cost to perturb
4406  double maximumFraction = 1.0e-5;
4407  int maxLength = 0;
4408  int minLength = numberRows_;
4409  double averageCost = 0.0;
4410  int numberNonZero = 0;
4411  // See if we need to perturb
4412  double *COIN_RESTRICT sort = new double[numberColumns_];
4413  // Use objective BEFORE scaling ??
4414  const double *COIN_RESTRICT obj = objective();
4415  for (int i = 0; i < numberColumns_; i++) {
4416    double value = fabs(obj[i]);
4417    sort[i] = value;
4418    averageCost += fabs(abcCost_[i + maximumAbcNumberRows_]);
4419    if (value)
4420      numberNonZero++;
4421  }
4422  if (numberNonZero)
4423    averageCost /= static_cast< double >(numberNonZero);
4424  else
4425    averageCost = 1.0;
4426  // allow for cost scaling
4427  averageCost *= objectiveScale_;
4428  std::sort(sort, sort + numberColumns_);
4429  int number = 1;
4430  double last = sort[0];
4431  for (int i = 1; i < numberColumns_; i++) {
4432    if (last != sort[i])
4433      number++;
4434    last = sort[i];
4435  }
4436  delete[] sort;
4437#ifdef PERT_STATISTICS
4438  printf("%d different values\n", number);
4439#endif
4440#if 0
4441  printf("nnz %d percent %d", number, (number * 100) / numberColumns_);
4442  if (number * 4 > numberColumns_)
4443    printf(" - Would not perturb\n");
4444  else
4445    printf(" - Would perturb\n");
4446#endif
4447  // If small costs then more dangerous (better to scale)
4448  double numberMultiplier = 1.25;
4449  if (averageCost < 1.0e-3)
4450    numberMultiplier = 20.0;
4451  else if (averageCost < 1.0e-2)
4452    numberMultiplier = 10.0;
4453  else if (averageCost < 1.0e-1)
4454    numberMultiplier = 3.0;
4455  else if (averageCost < 1.0)
4456    numberMultiplier = 2.0;
4457  numberMultiplier *= 2.0; // try
4458  if (number * numberMultiplier > numberColumns_ || perturbation_ >= 100) {
4459    perturbation_ = 101;
4460#if ABC_NORMAL_DEBUG > 0
4461    printf("Not perturbing - %d different costs, average %g\n",
4462      number, averageCost);
4463#endif
4464#ifdef CLEANUP_DJS
4465    // copy cost to perturbation
4466    CoinAbcMemcpy(abcPerturbation_, abcCost_, numberTotal_);
4467#endif
4468    return; // good enough
4469  }
4470  const int *COIN_RESTRICT columnLength = matrix_->getVectorLengths() - maximumAbcNumberRows_;
4471  const double *COIN_RESTRICT lower = lowerSaved_;
4472  const double *COIN_RESTRICT upper = upperSaved_;
4473  for (int iSequence = maximumAbcNumberRows_; iSequence < numberTotal_; iSequence++) {
4474    if (lower[iSequence] < upper[iSequence]) {
4475      int length = columnLength[iSequence];
4476      if (length > 2) {
4477        maxLength = CoinMax(maxLength, length);
4478        minLength = CoinMin(minLength, length);
4479      }
4480    }
4481  }
4482  bool uniformChange = false;
4483  bool inCbcOrOther = (specialOptions_ & 0x03000000) != 0;
4484  double constantPerturbation = 10.0 * dualTolerance_;
4485#ifdef PERT_STATISTICS
4486  breakdown("Objective before perturbation", numberColumns_, abcCost_ + numberRows_);
4487#endif
4488  //#define PERTURB_OLD_WAY
4489#ifndef PERTURB_OLD_WAY
4490#if 1
4491#ifdef HEAVY_PERTURBATION
4492  if (perturbation_ == 50)
4493    perturbation_ = HEAVY_PERTURBATION;
4494#endif
4495  if (perturbation_ >= 50) {
4496    // Experiment
4497    // maximumFraction could be 1.0e-10 to 1.0
4498    double m[] = { 1.0e-10, 1.0e-9, 1.0e-8, 1.0e-7, 1.0e-6, 1.0e-5, 1.0e-4, 1.0e-3, 1.0e-2, 1.0e-1, 1.0 };
4499    int whichOne = CoinMin(perturbation_ - 51, 10);
4500    if (whichOne < 0) {
4501      if (numberRows_ < 5000)
4502        whichOne = 2; // treat 50 as if 53
4503      else
4504        whichOne = 2; // treat 50 as if 52
4505    }
4506    if (inCbcOrOther && whichOne > 0)
4507      whichOne--;
4508    maximumFraction = m[whichOne];
4509    if (whichOne > 7)
4510      constantPerturbation *= 10.0;
4511    //if (whichOne<2)
4512    //constantPerturbation = maximumFraction*1.0e-10 * dualTolerance_;
4513  } else if (inCbcOrOther) {
4514    maximumFraction = 1.0e-6;
4515  }
4516#endif
4517  double smallestNonZero = 1.0e100;
4518  bool modifyRowCosts = false;
4519  numberNonZero = 0;
4520  //perturbation = 1.0e-10;
4521  perturbation = CoinMax(1.0e-10, maximumFraction);
4522  double *COIN_RESTRICT cost = abcCost_;
4523  bool allSame = true;
4524  double lastValue = 0.0;
4525
4526  for (int iRow = 0; iRow < numberRows_; iRow++) {
4527    double lo = lower[iRow];
4528    double up = upper[iRow];
4529    if (lo < up) {
4530      double value = fabs(cost[iRow]);
4531      perturbation = CoinMax(perturbation, value);
4532      if (value) {
4533        modifyRowCosts = true;
4534        smallestNonZero = CoinMin(smallestNonZero, value);
4535      }
4536    }
4537    if (lo && lo > -1.0e10) {
4538      numberNonZero++;
4539      lo = fabs(lo);
4540      if (!lastValue)
4541        lastValue = lo;
4542      else if (fabs(lo - lastValue) > 1.0e-7)
4543        allSame = false;
4544    }
4545    if (up && up < 1.0e10) {
4546      numberNonZero++;
4547      up = fabs(up);
4548      if (!lastValue)
4549        lastValue = up;
4550      else if (fabs(up - lastValue) > 1.0e-7)
4551        allSame = false;
4552    }
4553  }
4554  double lastValue2 = 0.0;
4555  for (int iSequence = maximumAbcNumberRows_; iSequence < numberTotal_; iSequence++) {
4556    double lo = lower[iSequence];
4557    double up = upper[iSequence];
4558    if (lo < up) {
4559      double value = fabs(cost[iSequence]);
4560      perturbation = CoinMax(perturbation, value);
4561      if (value) {
4562        smallestNonZero = CoinMin(smallestNonZero, value);
4563      }
4564    }
4565    if (lo && lo > -1.0e10) {
4566      //numberNonZero++;
4567      lo = fabs(lo);
4568      if (!lastValue2)
4569        lastValue2 = lo;
4570      else if (fabs(lo - lastValue2) > 1.0e-7)
4571        allSame = false;
4572    }
4573    if (up && up < 1.0e10) {
4574      //numberNonZero++;
4575      up = fabs(up);
4576      if (!lastValue2)
4577        lastValue2 = up;
4578      else if (fabs(up - lastValue2) > 1.0e-7)
4579        allSame = false;
4580    }
4581  }
4582  if (allSame) {
4583    // Check elements
4584    double smallestNegative;
4585    double largestNegative;
4586    double smallestPositive;
4587    double largestPositive;
4588    matrix_->rangeOfElements(smallestNegative, largestNegative,
4589      smallestPositive, largestPositive);
4590    if (smallestNegative == largestNegative && smallestPositive == largestPositive) {
4591      // Really hit perturbation
4592      double adjust = CoinMin(100.0 * maximumFraction, 1.0e-3 * CoinMax(lastValue, lastValue2));
4593      maximumFraction = CoinMax(adjust, maximumFraction);
4594    }
4595  }
4596  perturbation = CoinMin(perturbation, smallestNonZero / maximumFraction);
4597  double largestZero = 0.0;
4598  double largest = 0.0;
4599  double largestPerCent = 0.0;
4600  // modify costs
4601  bool printOut = (handler_->logLevel() == 63);
4602  printOut = false;
4603  //assert (!modifyRowCosts);
4604  modifyRowCosts = false;
4605  const double *COIN_RESTRICT perturbationArray = perturbationSaved_;
4606  if (modifyRowCosts) {
4607    for (int iRow = 0; iRow < numberRows_; iRow++) {
4608      if (lower[iRow] < upper[iRow]) {
4609        double value = perturbation;
4610        double currentValue = cost[iRow];
4611        value = CoinMin(value, maximumFraction * (fabs(currentValue) + 1.0e-1 * perturbation + 1.0e-3));
4612        double perturbationValue = 2.0 * (perturbationArray[iRow] - 0.5); // for now back to normal random
4613        if (lower[iRow] > -largeValue_) {
4614          if (fabs(lower[iRow]) < fabs(upper[iRow]))
4615            value *= perturbationValue;
4616          else
4617            value *= -perturbationValue;
4618        } else if (upper[iRow] < largeValue_) {
4619          value *= -perturbationValue;
4620        } else {
4621          value = 0.0;
4622        }
4623        if (currentValue) {
4624          largest = CoinMax(largest, fabs(value));
4625          if (fabs(value) > fabs(currentValue) * largestPerCent)
4626            largestPerCent = fabs(value / currentValue);
4627        } else {
4628          largestZero = CoinMax(largestZero, fabs(value));
4629        }
4630        if (printOut)
4631          printf("row %d cost %g change %g\n", iRow, cost[iRow], value);
4632        cost[iRow] += value;
4633      }
4634    }
4635  }
4636
4637  //double extraWeight = 10.0;
4638  // more its but faster double weight[]={1.0e-4,1.0e-2,1.0e-1,1.0,2.0,10.0,100.0,200.0,400.0,600.0,1000.0};
4639  // good its double weight[]={1.0e-4,1.0e-2,5.0e-1,1.0,2.0,5.0,10.0,20.0,30.0,40.0,100.0};
4640  double weight[] = { 1.0e-4, 1.0e-2, 5.0e-1, 1.0, 2.0, 5.0, 10.0, 20.0, 30.0, 40.0, 100.0 };
4641  //double weight[]={1.0e-4,1.0e-2,5.0e-1,1.0,20.0,50.0,100.0,120.0,130.0,140.0,200.0};
4642  // Scale back if wanted
4643  double weight2[] = { 1.0e-4, 1.0e-2, 5.0e-1, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0 };
4644  if (constantPerturbation < 99.0 * dualTolerance_ && false) {
4645    perturbation *= 0.1;
4646    //extraWeight = 0.5;
4647    memcpy(weight, weight2, sizeof(weight2));
4648  }
4649  // adjust weights if all columns long
4650  double factor = 1.0;
4651  if (maxLength) {
4652    factor = 3.0 / static_cast< double >(minLength);
4653  }
4654  // Make variables with more elements more expensive
4655  const double m1 = 0.5;
4656  double smallestAllowed = overallMultiplier * CoinMin(1.0e-2 * dualTolerance_, maximumFraction);
4657  double largestAllowed = overallMultiplier * CoinMax(1.0e3 * dualTolerance_, maximumFraction * averageCost);
4658  // smaller if in BAB
4659  //if (inCbcOrOther)
4660  //largestAllowed=CoinMin(largestAllowed,1.0e-5);
4661  //smallestAllowed = CoinMin(smallestAllowed,0.1*largestAllowed);
4662  randomNumberGenerator_.setSeed(1000000000);
4663#if ABC_NORMAL_DEBUG > 0
4664  printf("perturbation %g constantPerturbation %g smallestNonZero %g maximumFraction %g smallestAllowed %g\n",
4665    perturbation, constantPerturbation, smallestNonZero, maximumFraction,
4666    smallestAllowed);
4667#endif
4668  for (int iSequence = maximumAbcNumberRows_; iSequence < numberTotal_; iSequence++) {
4669    if (lower[iSequence] < upper[iSequence] && getInternalStatus(iSequence) != basic) {
4670      double value = perturbation;
4671      double currentValue = cost[iSequence];
4672      double currentLargestAllowed = fabs(currentValue) * 1.0e-3;
4673      if (!currentValue)
4674        currentLargestAllowed = 1.0e-3;
4675      if (currentLargestAllowed <= 10.0 * smallestAllowed)
4676        continue; // don't perturb tiny value
4677      value = CoinMin(value, constantPerturbation + maximumFraction * (fabs(currentValue) + 1.0e-1 * perturbation + 1.0e-8));
4678      //value = CoinMin(value,constantPerturbation;+maximumFraction*fabs(currentValue));
4679      double value2 = constantPerturbation + 1.0e-1 * smallestNonZero;
4680      if (uniformChange) {
4681        value = maximumFraction;
4682        value2 = maximumFraction;
4683      }
4684      //double perturbationValue=2.0*(perturbationArray[iSequence]-0.5); // for now back to normal random
4685      double perturbationValue = randomNumberGenerator_.randomDouble();
4686      perturbationValue = 1.0 - m1 + m1 * perturbationValue; // clean up
4687      if (lower[iSequence] > -largeValue_) {
4688        if (fabs(lower[iSequence]) < fabs(upper[iSequence])) {
4689          value *= perturbationValue;
4690          value2 *= perturbationValue;
4691        } else {
4692          //value *= -(1.0-m1+m1*randomNumberGenerator_.randomDouble());
4693          //value2 *= -(1.0-m1+m1*randomNumberGenerator_.randomDouble());
4694          value = 0.0;
4695        }
4696      } else if (upper[iSequence] < largeValue_) {
4697        value *= -perturbationValue;
4698        value2 *= -perturbationValue;
4699      } else {
4700        value = 0.0;
4701      }
4702      if (value) {
4703        int length = columnLength[iSequence];
4704        if (length > 3) {
4705          length = static_cast< int >(static_cast< double >(length) * factor);
4706          length = CoinMax(3, length);
4707        }
4708        double multiplier;
4709        if (length < 10)
4710          multiplier = weight[length];
4711        else
4712          multiplier = weight[10];
4713        value *= multiplier;
4714        value = CoinMin(value, value2);
4715        if (value) {
4716          // get in range
4717          if (fabs(value) <= smallestAllowed) {
4718            value *= 10.0;
4719            while (fabs(value) <= smallestAllowed)
4720              value *= 10.0;
4721          } else if (fabs(value) > currentLargestAllowed) {
4722            value *= 0.1;
4723            while (fabs(value) > currentLargestAllowed)
4724              value *= 0.1;
4725          }
4726        }
4727        if (currentValue) {
4728          largest = CoinMax(largest, fabs(value));
4729          if (fabs(value) > fabs(currentValue) * largestPerCent)
4730            largestPerCent = fabs(value / currentValue);
4731        } else {
4732          largestZero = CoinMax(largestZero, fabs(value));
4733        }
4734        // but negative if at ub
4735        if (getInternalStatus(iSequence) == atUpperBound)
4736          value = -value;
4737        if (printOut)
4738          printf("col %d cost %g change %g\n", iSequence, cost[iSequence], value);
4739        cost[iSequence] += value;
4740      }
4741    }
4742  }
4743#else
4744  // old way
4745  if (perturbation_ > 50) {
4746    // Experiment
4747    // maximumFraction could be 1.0e-10 to 1.0
4748    double m[] = { 1.0e-10, 1.0e-9, 1.0e-8, 1.0e-7, 1.0e-6, 1.0e-5, 1.0e-4, 1.0e-3, 1.0e-2, 1.0e-1, 1.0 };
4749    int whichOne = perturbation_ - 51;
4750    //if (inCbcOrOther&&whichOne>0)
4751    //whichOne--;
4752    maximumFraction = m[CoinMin(whichOne, 10)];
4753  } else if (inCbcOrOther) {
4754    //maximumFraction = 1.0e-6;
4755  }
4756  int iRow;
4757  double smallestNonZero = 1.0e100;
4758  numberNonZero = 0;
4759  if (perturbation_ >= 50) {
4760    perturbation = 1.0e-8;
4761    if (perturbation_ > 50 && perturbation_ < 60)
4762      perturbation = CoinMax(1.0e-8, maximumFraction);
4763    bool allSame = true;
4764    double lastValue = 0.0;
4765    for (iRow = 0; iRow < numberRows_; iRow++) {
4766      double lo = abcLower_[iRow];
4767      double up = abcUpper_[iRow];
4768      if (lo && lo > -1.0e10) {
4769        numberNonZero++;
4770        lo = fabs(lo);
4771        if (!lastValue)
4772          lastValue = lo;
4773        else if (fabs(lo - lastValue) > 1.0e-7)
4774          allSame = false;
4775      }
4776      if (up && up < 1.0e10) {
4777        numberNonZero++;
4778        up = fabs(up);
4779        if (!lastValue)
4780          lastValue = up;
4781        else if (fabs(up - lastValue) > 1.0e-7)
4782          allSame = false;
4783      }
4784    }
4785    double lastValue2 = 0.0;
4786    for (int iColumn = maximumAbcNumberRows_; iColumn < numberTotal_; iColumn++) {
4787      double lo = abcLower_[iColumn];
4788      double up = abcUpper_[iColumn];
4789      if (lo < up) {
4790        double value = fabs(abcCost_[iColumn]);
4791        perturbation = CoinMax(perturbation, value);
4792        if (value) {
4793          smallestNonZero = CoinMin(smallestNonZero, value);
4794        }
4795      }
4796      if (lo && lo > -1.0e10) {
4797        //numberNonZero++;
4798        lo = fabs(lo);
4799        if (!lastValue2)
4800          lastValue2 = lo;
4801        else if (fabs(lo - lastValue2) > 1.0e-7)
4802          allSame = false;
4803      }
4804      if (up && up < 1.0e10) {
4805        //numberNonZero++;
4806        up = fabs(up);
4807        if (!lastValue2)
4808          lastValue2 = up;
4809        else if (fabs(up - lastValue2) > 1.0e-7)
4810          allSame = false;
4811      }
4812    }
4813    if (allSame) {
4814      // Check elements
4815      double smallestNegative;
4816      double largestNegative;
4817      double smallestPositive;
4818      double largestPositive;
4819      matrix_->rangeOfElements(smallestNegative, largestNegative,
4820        smallestPositive, largestPositive);
4821      if (smallestNegative == largestNegative && smallestPositive == largestPositive) {
4822        // Really hit perturbation
4823        double adjust = CoinMin(100.0 * maximumFraction, 1.0e-3 * CoinMax(lastValue, lastValue2));
4824        maximumFraction = CoinMax(adjust, maximumFraction);
4825      }
4826    }
4827    perturbation = CoinMin(perturbation, smallestNonZero / maximumFraction);
4828  }
4829  double largestZero = 0.0;
4830  double largest = 0.0;
4831  double largestPerCent = 0.0;
4832  // modify costs
4833  bool printOut = (handler_->logLevel() == 63);
4834  printOut = false;
4835  // more its but faster double weight[]={1.0e-4,1.0e-2,1.0e-1,1.0,2.0,10.0,100.0,200.0,400.0,600.0,1000.0};
4836  // good its double weight[]={1.0e-4,1.0e-2,5.0e-1,1.0,2.0,5.0,10.0,20.0,30.0,40.0,100.0};
4837  double weight[] = { 1.0e-4, 1.0e-2, 5.0e-1, 1.0, 2.0, 5.0, 10.0, 20.0, 30.0, 40.0, 100.0 };
4838  //double weight[]={1.0e-4,1.0e-2,5.0e-1,1.0,20.0,50.0,100.0,120.0,130.0,140.0,200.0};
4839  //double extraWeight = 10.0;
4840  // Scale back if wanted
4841  double weight2[] = { 1.0e-4, 1.0e-2, 5.0e-1, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0 };
4842  if (constantPerturbation < 99.0 * dualTolerance_) {
4843    perturbation *= 0.1;
4844    //extraWeight = 0.5;
4845    memcpy(weight, weight2, sizeof(weight2));
4846  }
4847  // adjust weights if all columns long
4848  double factor = 1.0;
4849  if (maxLength) {
4850    factor = 3.0 / static_cast< double >(minLength);
4851  }
4852  // Make variables with more elements more expensive
4853  const double m1 = 0.5;
4854  double smallestAllowed = CoinMin(1.0e-2 * dualTolerance_, maximumFraction);
4855  double largestAllowed = CoinMax(1.0e3 * dualTolerance_, maximumFraction * averageCost);
4856  // smaller if in BAB
4857  //if (inCbcOrOther)
4858  //largestAllowed=CoinMin(largestAllowed,1.0e-5);
4859  //smallestAllowed = CoinMin(smallestAllowed,0.1*largestAllowed);
4860  for (int iColumn = maximumAbcNumberRows_; iColumn < numberTotal_; iColumn++) {
4861    if (abcLower_[iColumn] < abcUpper_[iColumn]
4862      && getInternalStatus(iColumn) != basic) {
4863      double value = perturbation;
4864      double currentValue = abcCost_[iColumn];
4865      value = CoinMin(value, constantPerturbation + maximumFraction * (fabs(currentValue) + 1.0e-1 * perturbation + 1.0e-8));
4866      //value = CoinMin(value,constantPerturbation;+maximumFraction*fabs(currentValue));
4867      double value2 = constantPerturbation + 1.0e-1 * smallestNonZero;
4868      if (uniformChange) {
4869        value = maximumFraction;
4870        value2 = maximumFraction;
4871      }
4872      if (abcLower_[iColumn] > -largeValue_) {
4873        if (fabs(abcLower_[iColumn]) < fabs(abcUpper_[iColumn])) {
4874          value *= (1.0 - m1 + m1 * randomNumberGenerator_.randomDouble());
4875          value2 *= (1.0 - m1 + m1 * randomNumberGenerator_.randomDouble());
4876        } else {
4877          //value *= -(1.0-m1+m1*randomNumberGenerator_.randomDouble());
4878          //value2 *= -(1.0-m1+m1*randomNumberGenerator_.randomDouble());
4879          value = 0.0;
4880        }
4881      } else if (abcUpper_[iColumn] < largeValue_) {
4882        value *= -(1.0 - m1 + m1 * randomNumberGenerator_.randomDouble());
4883        value2 *= -(1.0 - m1 + m1 * randomNumberGenerator_.randomDouble());
4884      } else {
4885        value = 0.0;
4886      }
4887      if (value) {
4888        int length = columnLength[iColumn];
4889        if (length > 3) {
4890          length = static_cast< int >(static_cast< double >(length) * factor);
4891          length = CoinMax(3, length);
4892        }
4893        double multiplier;
4894#if 1
4895        if (length < 10)
4896          multiplier = weight[length];
4897        else
4898          multiplier = weight[10];
4899#else
4900        if (length < 10)
4901          multiplier = weight[length];
4902        else
4903          multiplier = weight[10] + extraWeight * (length - 10);
4904        multiplier *= 0.5;
4905#endif
4906        value *= multiplier;
4907        value = CoinMin(value, value2);
4908        if (value) {
4909          // get in range
4910          if (fabs(value) <= smallestAllowed) {
4911            value *= 10.0;
4912            while (fabs(value) <= smallestAllowed)
4913              value *= 10.0;
4914          } else if (fabs(value) > largestAllowed) {
4915            value *= 0.1;
4916            while (fabs(value) > largestAllowed)
4917              value *= 0.1;
4918          }
4919        }
4920        if (currentValue) {
4921          largest = CoinMax(largest, fabs(value));
4922          if (fabs(value) > fabs(currentValue) * largestPerCent)
4923            largestPerCent = fabs(value / currentValue);
4924        } else {
4925          largestZero = CoinMax(largestZero, fabs(value));
4926        }
4927        // but negative if at ub
4928        if (getInternalStatus(iColumn) == atUpperBound)
4929          value = -value;
4930        if (printOut)
4931          printf("col %d cost %g change %g\n", iColumn, objectiveWork_[iColumn], value);
4932        abcCost_[iColumn] += value;
4933      }
4934    }
4935  }
4936#endif
4937#ifdef PERT_STATISTICS
4938  {
4939    double averageCost = 0.0;
4940    int numberNonZero = 0;
4941    double *COIN_RESTRICT sort = new double[numberColumns_];
4942    for (int i = 0; i < numberColumns_; i++) {
4943      double value = fabs(abcCost_[i + maximumAbcNumberRows_]);
4944      sort[i] = value;
4945      averageCost += value;
4946      if (value)
4947        numberNonZero++;
4948    }
4949    if (numberNonZero)
4950      averageCost /= static_cast< double >(numberNonZero);
4951    else
4952      averageCost = 1.0;
4953    std::sort(sort, sort + numberColumns_);
4954    int number = 1;
4955    double last = sort[0];
4956    for (int i = 1; i < numberColumns_; i++) {
4957      if (last != sort[i])
4958        number++;
4959      last = sort[i];
4960    }
4961    delete[] sort;
4962#if ABC_NORMAL_DEBUG > 0
4963    printf("nnz %d percent %d", number, (number * 100) / numberColumns_);
4964    breakdown("Objective after perturbation", numberColumns_, abcCost_ + numberRows_);
4965#endif
4966  }
4967#endif
4968  perturbation_ = 100;
4969  // copy cost to  perturbation
4970#ifdef CLEANUP_DJS
4971  CoinAbcMemcpy(abcPerturbation_, abcCost_, numberTotal_);
4972#endif
4973  handler_->message(CLP_SIMPLEX_PERTURB, messages_)
4974    << 100.0 * maximumFraction << perturbation << largest << 100.0 * largestPerCent << largestZero
4975    << CoinMessageEol;
4976}
4977/* Does something about fake tolerances
4978   -1 initializes
4979   1 might be optimal (back to original costs)
4980   2 after change currentDualBound
4981   3 might be optimal BUT change currentDualBound
4982*/
4983int AbcSimplexDual::bounceTolerances(int type)
4984{
4985  // for now turn off
4986  if (type == -1) {
4987    if (!perturbationFactor_)
4988      perturbationFactor_ = 0.001;
4989    CoinAbcMemcpy(abcCost_, costSaved_, numberTotal_);
4990    if (!perturbationFactor_) {
4991      //currentDualBound_=CoinMin(dualBound_,1.0e11);
4992      //CoinAbcMultiplyAdd(perturbationSaved_,numberTotal_,currentDualTolerance_,
4993      //                 abcPerturbation_,0.0);
4994    } else {
4995      //currentDualBound_=CoinMin(dualBound_,1.0e11);
4996      double factor = currentDualTolerance_ * perturbationFactor_;
4997      perturbB(factor, 0);
4998    }
4999#if ABC_NORMAL_DEBUG > 0
5000    //printf("new dual bound %g (bounceTolerances)\n",currentDualBound_);
5001#endif
5002  } else {
5003    if (stateDualColumn_ == -2) {
5004      double newTolerance = CoinMin(1.0e-1, dualTolerance_ * 100.0);
5005#if ABC_NORMAL_DEBUG > 0
5006      printf("After %d iterations reducing dual tolerance from %g to %g\n",
5007        numberIterations_, currentDualTolerance_, newTolerance);
5008#endif
5009      currentDualTolerance_ = newTolerance;
5010      stateDualColumn_ = -1;
5011    } else if (stateDualColumn_ == -1) {
5012#if ABC_NORMAL_DEBUG > 0
5013      printf("changing dual tolerance from %g to %g (vaguely optimal)\n",
5014        currentDualTolerance_, dualTolerance_);
5015#endif
5016      currentDualTolerance_ = dualTolerance_;
5017      stateDualColumn_ = 0;
5018    } else {
5019      if (stateDualColumn_ > 1) {
5020        normalDualColumnIteration_ = numberIterations_ + CoinMax(10, numberRows_ / 4);
5021#if ABC_NORMAL_DEBUG > 0
5022        printf("no cost modification until iteration %d\n", normalDualColumnIteration_);
5023#endif
5024      }
5025    }
5026    resetFakeBounds(0);
5027    CoinAbcMemcpy(abcCost_, abcPerturbation_, numberTotal_);
5028    numberChanged_ = 0; // Number of variables with changed costs
5029    moveToBasic();
5030    if (type != 1)
5031      return 0;
5032    // save status and solution
5033    CoinAbcMemcpy(status_, internalStatus_, numberTotal_);
5034    CoinAbcMemcpy(solutionSaved_, abcSolution_, numberTotal_);
5035#if ABC_NORMAL_DEBUG > 0
5036    printf("A sumDual %g sumPrimal %g\n", sumDualInfeasibilities_, sumPrimalInfeasibilities_);
5037#endif
5038    gutsOfSolution(3);
5039#if ABC_NORMAL_DEBUG > 0
5040    printf("B sumDual %g sumPrimal %g\n", sumDualInfeasibilities_, sumPrimalInfeasibilities_);
5041#endif
5042    double saveSum1 = sumPrimalInfeasibilities_;
5043    double saveSumD1 = sumDualInfeasibilities_;
5044    bool goToPrimal = false;
5045    if (!sumOfRelaxedDualInfeasibilities_) {
5046      if (perturbation_ < 101) {
5047        // try going back to original costs
5048        CoinAbcMemcpy(abcCost_, costSaved_, numberTotal_);
5049        perturbation_ = 101;
5050        numberChanged_ = 0; // Number of variables with changed costs
5051        moveToBasic();
5052        gutsOfSolution(3);
5053#if ABC_NORMAL_DEBUG > 0
5054        printf("B2 sumDual %g sumPrimal %g\n", sumDualInfeasibilities_, sumPrimalInfeasibilities_);
5055#endif
5056        saveSumD1 = sumDualInfeasibilities_;
5057      }
5058      if (sumOfRelaxedDualInfeasibilities_) {
5059        makeNonFreeVariablesDualFeasible();
5060        moveToBasic();
5061        abcProgress_.reset();
5062        gutsOfSolution(3);
5063#if ABC_NORMAL_DEBUG > 0
5064        printf("C sumDual %g sumPrimal %g\n", sumDualInfeasibilities_, sumPrimalInfeasibilities_);
5065#endif
5066        if (sumPrimalInfeasibilities_ > 1.0e2 * saveSum1 + currentDualBound_ * 0.2
5067          /*|| sumOfRelaxedDualInfeasibilities_<1.0e-4*sumPrimalInfeasibilities_*/) {
5068          numberDisasters_++;
5069          if (numberDisasters_ > 2)
5070            goToPrimal = true;
5071        }
5072      } else {
5073        numberDualInfeasibilities_ = 0;
5074        sumDualInfeasibilities_ = 0.0;
5075      }
5076    } else {
5077      // carry on
5078      makeNonFreeVariablesDualFeasible();
5079      moveToBasic();
5080      abcProgress_.reset();
5081      gutsOfSolution(3);
5082#if ABC_NORMAL_DEBUG > 0
5083      printf("C2 sumDual %g sumPrimal %g\n", sumDualInfeasibilities_, sumPrimalInfeasibilities_);
5084#endif
5085      numberDisasters_++;
5086      if (sumPrimalInfeasibilities_ > 1.0e2 * saveSum1 + currentDualBound_ * 0.2
5087        || (saveSumD1 < 1.0e-4 * sumPrimalInfeasibilities_ && saveSumD1 < 1.0)) {
5088        numberDisasters_++;
5089        if (numberDisasters_ > 2 || saveSumD1 < 1.0e-10 * sumPrimalInfeasibilities_)
5090          goToPrimal = true;
5091      }
5092    }
5093    if (goToPrimal) {
5094      // go to primal
5095      CoinAbcMemcpy(internalStatus_, status_, numberTotal_);
5096      CoinAbcMemcpy(abcSolution_, solutionSaved_, numberTotal_);
5097      int *COIN_RESTRICT abcPivotVariable = abcPivotVariable_;
5098      // redo pivotVariable
5099      int numberBasic = 0;
5100      for (int i = 0; i < numberTotal_; i++) {
5101        if (getInternalStatus(i) == basic)
5102          abcPivotVariable[numberBasic++] = i;
5103      }
5104      assert(numberBasic == numberRows_);
5105      moveToBasic();
5106      checkPrimalSolution(false);
5107#if ABC_NORMAL_DEBUG > 0
5108      printf("D sumPrimal %g\n", sumPrimalInfeasibilities_);
5109#endif
5110      return 1;
5111    }
5112  }
5113
5114  CoinAbcGatherFrom(abcCost_, costBasic_, abcPivotVariable_, numberRows_);
5115  return 0;
5116}
5117#if ABC_NORMAL_DEBUG > 0
5118extern int cilk_conflict;
5119#endif
5120int AbcSimplexDual::noPivotRow()
5121{
5122#if ABC_NORMAL_DEBUG > 3
5123  if (handler_->logLevel() & 32)
5124    printf("** no row pivot\n");
5125#endif
5126  int numberPivots = abcFactorization_->pivots();
5127  bool specialCase;
5128  int useNumberFake;
5129  int returnCode = 0;
5130  if (numberPivots <= CoinMax(dontFactorizePivots_, 20) && (specialOptions_ & 2048) != 0 && currentDualBound_ >= 1.0e8) {
5131    specialCase = true;
5132    // as dual bound high - should be okay
5133    useNumberFake = 0;
5134  } else {
5135    specialCase = false;
5136    useNumberFake = numberAtFakeBound();
5137  }
5138  if (!numberPivots || specialCase) {
5139    // may have crept through - so may be optimal
5140    // check any flagged variables
5141    if (numberFlagged_ && numberPivots) {
5142      // try factorization
5143      returnCode = -2;
5144    }
5145
5146    if (useNumberFake || numberDualInfeasibilities_) {
5147      // may be dual infeasible
5148      if ((specialOptions_ & 1024) == 0)
5149        problemStatus_ = -5;
5150      else if (!useNumberFake && numberPrimalInfeasibilities_
5151        && !numberPivots)
5152        problemStatus_ = 1;
5153    } else {
5154      if (numberFlagged_) {
5155#if ABC_NORMAL_DEBUG > 3
5156        std::cout << "Flagged variables at end - infeasible?" << std::endl;
5157        printf("Probably infeasible - pivot was %g\n", alpha_);
5158#endif
5159#if ABC_NORMAL_DEBUG > 3
5160        abort();
5161#endif
5162        problemStatus_ = -5;
5163      } else {
5164        problemStatus_ = 0;
5165#ifndef CLP_CHECK_NUMBER_PIVOTS
5166#define CLP_CHECK_NUMBER_PIVOTS 10
5167#endif
5168#if CLP_CHECK_NUMBER_PIVOTS < 20
5169        if (numberPivots > CLP_CHECK_NUMBER_PIVOTS) {
5170          gutsOfSolution(2);
5171          rawObjectiveValue_ += sumNonBasicCosts_;
5172          setClpSimplexObjectiveValue();
5173          if (numberPrimalInfeasibilities_) {
5174            problemStatus_ = -1;
5175            returnCode = -2;
5176          }
5177        }
5178#endif
5179        if (!problemStatus_) {
5180          // make it look OK
5181          numberPrimalInfeasibilities_ = 0;
5182          sumPrimalInfeasibilities_ = 0.0;
5183          numberDualInfeasibilities_ = 0;
5184          sumDualInfeasibilities_ = 0.0;
5185          if (numberChanged_) {
5186            numberChanged_ = 0; // Number of variables with changed costs
5187            CoinAbcMemcpy(abcCost_, abcCost_ + maximumNumberTotal_, numberTotal_);
5188            CoinAbcGatherFrom(abcCost_, costBasic_, abcPivotVariable_, numberRows_);
5189            // make sure duals are current
5190            gutsOfSolution(1);
5191            setClpSimplexObjectiveValue();
5192            abcProgress_.modifyObjective(-COIN_DBL_MAX);
5193            if (numberDualInfeasibilities_) {
5194              problemStatus_ = 10; // was -3;
5195              numberTimesOptimal_++;
5196            } else {
5197              computeObjectiveValue(true);
5198            }
5199          } else if (numberPivots) {
5200            computeObjectiveValue(true);
5201          }
5202          sumPrimalInfeasibilities_ = 0.0;
5203        }
5204        if ((specialOptions_ & (1024 + 16384)) != 0 && !problemStatus_ && abcFactorization_->pivots()) {
5205          CoinIndexedVector *arrayVector = &usefulArray_[arrayForFtran_];
5206          double *rhs = arrayVector->denseVector();
5207          arrayVector->clear();
5208          CoinAbcScatterTo(solutionBasic_, abcSolution_, abcPivotVariable_, numberRows_);
5209          abcMatrix_->timesModifyExcludingSlacks(-1.0, abcSolution_, rhs);
5210          //#define CHECK_ACCURACY
5211#ifdef CHECK_ACCURACY
5212          bool bad = false;
5213#endif
5214          bool bad2 = false;
5215          int i;
5216          for (i = 0; i < numberRows_; i++) {
5217            if (rhs[i] < abcLower_[i] - primalTolerance_ || rhs[i] > abcUpper_[i] + primalTolerance_) {
5218              bad2 = true;
5219#ifdef CHECK_ACCURACY
5220              printf("row %d out of bounds %g, %g correct %g bad %g\n", i,
5221                abcLower_[i], abcUpper_[i],
5222                rhs[i], abcSolution_[i]);
5223#endif
5224            } else if (fabs(rhs[i] - abcSolution_[i]) > 1.0e-3) {
5225#ifdef CHECK_ACCURACY
5226              bad = true;
5227              printf("row %d correct %g bad %g\n", i, rhs[i], abcSolution_[i]);
5228#endif
5229            }
5230            rhs[i] = 0.0;
5231          }
5232          for (int i = 0; i < numberTotal_; i++) {
5233            if (abcSolution_[i] < abcLower_[i] - primalTolerance_ || abcSolution_[i] > abcUpper_[i] + primalTolerance_) {
5234              bad2 = true;
5235#ifdef CHECK_ACCURACY
5236              printf("column %d out of bounds %g, %g value %g\n", i,
5237                abcLower_[i], abcUpper_[i],
5238                abcSolution_[i]);
5239#endif
5240            }
5241          }
5242          if (bad2) {
5243            problemStatus_ = -3;
5244            returnCode = -2;
5245            // Force to re-factorize early next time
5246            int numberPivots = abcFactorization_->pivots();
5247            if (forceFactorization_ < 0)
5248              forceFactorization_ = 100000;
5249            forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
5250          }
5251        }
5252      }
5253    }
5254  } else {
5255    problemStatus_ = -3;
5256    returnCode = -2;
5257    // Force to re-factorize early next time
5258    int numberPivots = abcFactorization_->pivots();
5259    if (forceFactorization_ < 0)
5260      forceFactorization_ = 100000;
5261    forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
5262  }
5263  return returnCode;
5264}
5265void AbcSimplexDual::dualPivotColumn()
5266{
5267  double saveAcceptable = acceptablePivot_;
5268  if (largestPrimalError_ > 1.0e-5 || largestDualError_ > 1.0e-5) {
5269    //if ((largestPrimalError_>1.0e-5||largestDualError_>1.0e-5)&&false) {
5270    if (!abcFactorization_->pivots())
5271      acceptablePivot_ *= 1.0e2;
5272    else if (abcFactorization_->pivots() < 5)
5273      acceptablePivot_ *= 1.0e1;
5274  }
5275  dualColumn2();
5276  if (sequenceIn_ < 0 && acceptablePivot_ > saveAcceptable) {
5277    acceptablePivot_ = saveAcceptable;
5278#if ABC_NORMAL_DEBUG > 0
5279    printf("No pivot!!!\n");
5280#endif
5281    // try again
5282    dualColumn1(true);
5283    dualColumn2();
5284  } else {
5285    acceptablePivot_ = saveAcceptable;
5286  }
5287  if ((stateOfProblem_ & VALUES_PASS) != 0) {
5288    // check no flips
5289    assert(!usefulArray_[arrayForFlipBounds_].getNumElements());
5290    // If no pivot - then try other way if plausible
5291    if (sequenceIn_ < 0) {
5292      // could use fake primal basic values
5293      if ((directionOut_ < 0 && lowerOut_ > -1.0e30) || (directionOut_ > 0 && upperOut_ < 1.30)) {
5294        directionOut_ = -directionOut_;
5295        if (directionOut_ < 0 && abcDj_[sequenceOut_] < 0.0)
5296          upperTheta_ = -abcDj_[sequenceOut_];
5297        else if (directionOut_ > 0 && abcDj_[sequenceOut_] > 0.0)
5298          upperTheta_ = abcDj_[sequenceOut_];
5299        else
5300          upperTheta_ = 1.0e30;
5301        CoinPartitionedVector &tableauRow = usefulArray_[arrayForTableauRow_];
5302        CoinPartitionedVector &candidateList = usefulArray_[arrayForDualColumn_];
5303        assert(!candidateList.getNumElements());
5304        tableauRow.compact();
5305        candidateList.compact();
5306        // redo candidate list
5307        int *COIN_RESTRICT index = tableauRow.getIndices();
5308        double *COIN_RESTRICT array = tableauRow.denseVector();
5309        double *COIN_RESTRICT arrayCandidate = candidateList.denseVector();
5310        int *COIN_RESTRICT indexCandidate = candidateList.getIndices();
5311        const double *COIN_RESTRICT abcDj = abcDj_;
5312        const unsigned char *COIN_RESTRICT internalStatus = internalStatus_;
5313        // do first pass to get possibles
5314        double bestPossible = 0.0;
5315        // We can also see if infeasible or pivoting on free
5316        double tentativeTheta = 1.0e25; // try with this much smaller as guess
5317        double acceptablePivot = currentAcceptablePivot_;
5318        double dualT = -currentDualTolerance_;
5319        // fixed will have been taken out by now
5320        const double multiplier[] = { 1.0, -1.0 };
5321        int freeSequence = -1;
5322        double freeAlpha = 0.0;
5323        int numberNonZero = tableauRow.getNumElements();
5324        int numberRemaining = 0;
5325        if (ordinaryVariables()) {
5326          for (int i = 0; i < numberNonZero; i++) {
5327            int iSequence = index[i];
5328            double tableauValue = array[i];
5329            unsigned char iStatus = internalStatus[iSequence] & 7;
5330            double mult = multiplier[iStatus];
5331            double alpha = tableauValue * mult;
5332            double oldValue = abcDj[iSequence] * mult;
5333            double value = oldValue - tentativeTheta * alpha;
5334            if (value < dualT) {
5335              bestPossible = CoinMax(bestPossible, alpha);
5336              value = oldValue - upperTheta_ * alpha;
5337              if (value < dualT && alpha >= acceptablePivot) {
5338                upperTheta_ = (oldValue - dualT) / alpha;
5339              }
5340              // add to list
5341              arrayCandidate[numberRemaining] = alpha;
5342              indexCandidate[numberRemaining++] = iSequence;
5343            }
5344          }
5345        } else {
5346          double badFree = 0.0;
5347          freeAlpha = currentAcceptablePivot_;
5348          double currentDualTolerance = currentDualTolerance_;
5349          for (int i = 0; i < numberNonZero; i++) {
5350            int iSequence = index[i];
5351            double tableauValue = array[i];
5352            unsigned char iStatus = internalStatus[iSequence] & 7;
5353            if ((iStatus & 4) == 0) {
5354              double mult = multiplier[iStatus];
5355              double alpha = tableauValue * mult;
5356              double oldValue = abcDj[iSequence] * mult;
5357              double value = oldValue - tentativeTheta * alpha;
5358              if (value < dualT) {
5359                bestPossible = CoinMax(bestPossible, alpha);
5360                value = oldValue - upperTheta_ * alpha;
5361                if (value < dualT && alpha >= acceptablePivot) {
5362                  upperTheta_ = (oldValue - dualT) / alpha;
5363                }
5364                // add to list
5365                arrayCandidate[numberRemaining] = alpha;
5366                indexCandidate[numberRemaining++] = iSequence;
5367              }
5368            } else {
5369              bool keep;
5370              bestPossible = CoinMax(bestPossible, fabs(tableauValue));
5371              double oldValue = abcDj[iSequence];
5372              // If free has to be very large - should come in via dualRow
5373              //if (getInternalStatus(iSequence+addSequence)==isFree&&fabs(tableauValue)<1.0e-3)
5374              //break;
5375              if (oldValue > currentDualTolerance) {
5376                keep = true;
5377              } else if (oldValue < -currentDualTolerance) {
5378                keep = true;
5379              } else {
5380                if (fabs(tableauValue) > CoinMax(10.0 * acceptablePivot, 1.0e-5)) {
5381                  keep = true;
5382                } else {
5383                  keep = false;
5384                  badFree = CoinMax(badFree, fabs(tableauValue));
5385                }
5386              }
5387              if (keep) {
5388                // free - choose largest
5389                if (fabs(tableauValue) > fabs(freeAlpha)) {
5390                  freeAlpha = tableauValue;
5391                  freeSequence = iSequence;
5392                }
5393              }
5394            }
5395          }
5396        }
5397        candidateList.setNumElements(numberRemaining);
5398        if (freeSequence >= 0) {
5399          sequenceIn_ = freeSequence;
5400        } else if (fabs(upperTheta_ - fabs(abcDj_[sequenceOut_])) < dualTolerance_) {
5401          // can just move
5402          stateOfIteration_ = -1;
5403          return;
5404        }
5405        dualColumn2();
5406      }
5407    }
5408    // redo dualOut_
5409    if (directionOut_ < 0) {
5410      dualOut_ = valueOut_ - upperOut_;
5411    } else {
5412      dualOut_ = lowerOut_ - valueOut_;
5413    }
5414  }
5415  usefulArray_[arrayForDualColumn_].clear();
5416  //clearArrays(arrayForDualColumn_);
5417  if (sequenceIn_ >= 0 && upperTheta_ == COIN_DBL_MAX) {
5418    lowerIn_ = abcLower_[sequenceIn_];
5419    upperIn_ = abcUpper_[sequenceIn_];
5420    valueIn_ = abcSolution_[sequenceIn_];
5421    dualIn_ = abcDj_[sequenceIn_];
5422    if (valueIn_ < lowerIn_ + primalTolerance_) {
5423      directionIn_ = 1;
5424    } else if (valueIn_ > upperIn_ - primalTolerance_) {
5425      directionIn_ = -1;
5426    } else {
5427      if (alpha_ < 0.0) {
5428        // as if from upper bound
5429        directionIn_ = -1;
5430        //assert(upperIn_ == valueIn_);
5431      } else {
5432        // as if from lower bound
5433        directionIn_ = 1;
5434        //assert(lowerIn_ == valueIn_);
5435      }
5436    }
5437  }
5438#if ABC_DEBUG
5439  checkArrays(4);
5440#endif
5441#if CAN_HAVE_ZERO_OBJ > 1
5442  if ((specialOptions_ & 2097152) != 0)
5443    theta_ = 0.0;
5444#endif
5445  movement_ = 0.0;
5446  objectiveChange_ = 0.0;
5447  /*
5448    0 - take iteration
5449    1 - don't take but continue
5450    2 - don't take and break
5451  */
5452  btranAlpha_ = -alpha_ * directionOut_; // for check
5453}
5454void AbcSimplexDual::getTableauColumnPart2()
5455{
5456#if ABC_PARALLEL
5457  if (parallelMode_ != 0) {
5458    abcFactorization_->updateColumnFTPart2(usefulArray_[arrayForFtran_]);
5459    // pivot element
5460    alpha_ = usefulArray_[arrayForFtran_].denseVector()[pivotRow_];
5461  }
5462#endif
5463}
5464void AbcSimplexDual::replaceColumnPart3()
5465{
5466  abcFactorization_->replaceColumnPart3(this,
5467    &usefulArray_[arrayForReplaceColumn_],
5468    &usefulArray_[arrayForFtran_],
5469    lastPivotRow_,
5470    ftAlpha_ ? ftAlpha_ : alpha_);
5471}
5472void AbcSimplexDual::checkReplacePart1()
5473{
5474  //abcFactorization_->checkReplacePart1a(&usefulArray_[arrayForReplaceColumn_],pivotRow_);
5475  //usefulArray_[arrayForReplaceColumn_].print();
5476  ftAlpha_ = abcFactorization_->checkReplacePart1(&usefulArray_[arrayForReplaceColumn_], pivotRow_);
5477}
5478#if 0
5479void
5480AbcSimplexDual::checkReplacePart1a()
5481{
5482  abcFactorization_->checkReplacePart1a(&usefulArray_[arrayForReplaceColumn_],pivotRow_);
5483}
5484void
5485AbcSimplexDual::checkReplacePart1b()
5486{
5487  //usefulArray_[arrayForReplaceColumn_].print();
5488  ftAlpha_=abcFactorization_->checkReplacePart1b(&usefulArray_[arrayForReplaceColumn_],pivotRow_);
5489}
5490#endif
5491/*
5492  returns
5493  0 - OK
5494  1 - take iteration then re-factorize
5495  2 - flag something and skip
5496  3 - break and re-factorize
5497  5 - take iteration then re-factorize because of memory
5498 */
5499int AbcSimplexDual::checkReplace()
5500{
5501  int returnCode = 0;
5502  if (!stateOfIteration_) {
5503    int saveStatus = problemStatus_;
5504    // check update
5505#if MOVE_REPLACE_PART1A < 0
5506    checkReplacePart1();
5507#endif
5508    int updateStatus = abcFactorization_->checkReplacePart2(pivotRow_, btranAlpha_, alpha_, ftAlpha_);
5509#if ABC_DEBUG
5510    checkArrays();
5511#endif
5512    // If looks like bad pivot - refactorize
5513    if (fabs(dualOut_) > 1.0e50)
5514      updateStatus = 2;
5515    // if no pivots, bad update but reasonable alpha - take and invert
5516    if (updateStatus == 2 && !abcFactorization_->pivots() && fabs(alpha_) > 1.0e-5)
5517      updateStatus = 4;
5518    if (updateStatus == 1 || updateStatus == 4) {
5519      // slight error
5520      if (abcFactorization_->pivots() > 5 || updateStatus == 4) {
5521        problemStatus_ = -2; // factorize now
5522        returnCode = -3;
5523      }
5524    } else if (updateStatus == 2) {
5525      // major error
5526</