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

Last change on this file since 2030 was 2024, checked in by forrest, 5 years ago

changes to abc

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