source: stable/1.15/Clp/src/AbcSimplexDual.cpp @ 1989

Last change on this file since 1989 was 1989, checked in by forrest, 6 years ago

Minor stability fixes and an option to allow perturbation after presolve

  • Property svn:keywords set to Id
File size: 201.9 KB
<
Line 
1/* $Id: AbcSimplexDual.cpp 1989 2013-11-11 17:27:32Z 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  if (primalFeasible()&&ignoreStuff!=2) {
3297    // may be optimal - or may be bounds are wrong
3298    int numberChangedBounds = changeBounds(0, changeCost);
3299    if (numberChangedBounds <= 0 && !numberDualInfeasibilities_) {
3300      //looks optimal - do we need to reset tolerance
3301      if (lastCleaned_ < numberIterations_ && /*numberTimesOptimal_*/ stateDualColumn_ < 4 &&
3302          (numberChanged_ || (specialOptions_ & 4096) == 0)) {
3303#if CLP_CAN_HAVE_ZERO_OBJ
3304        if ((specialOptions_&2097152)==0) {
3305#endif
3306          numberTimesOptimal_++;
3307          if (stateDualColumn_==-2) {
3308#if ABC_NORMAL_DEBUG>0
3309            printf("Going straight from state -2 to 0\n");
3310#endif
3311            stateDualColumn_=-1;
3312          }
3313          changeMade_++; // say something changed
3314          if (numberTimesOptimal_ == 1) {
3315            if (fabs(currentDualTolerance_-dualTolerance_)>1.0e-12) {
3316#if ABC_NORMAL_DEBUG>0
3317              printf("changing dual tolerance from %g to %g (numberTimesOptimal_==%d)\n",
3318                     currentDualTolerance_,dualTolerance_,numberTimesOptimal_);
3319#endif
3320              currentDualTolerance_ = dualTolerance_;
3321            }
3322          } else {
3323            if (numberTimesOptimal_ == 2) {
3324              // better to have small tolerance even if slower
3325              abcFactorization_->zeroTolerance(CoinMin(abcFactorization_->zeroTolerance(), 1.0e-15));
3326            }
3327#if ABC_NORMAL_DEBUG>0
3328            printf("changing dual tolerance from %g to %g (numberTimesOptimal_==%d)\n",
3329                   currentDualTolerance_,dualTolerance_*pow(2.0,numberTimesOptimal_-1),numberTimesOptimal_);
3330#endif
3331            currentDualTolerance_ = dualTolerance_ * pow(2.0, numberTimesOptimal_ - 1);
3332          }
3333          if (numberChanged_>0) {
3334            handler_->message(CLP_DUAL_ORIGINAL, messages_)
3335              << CoinMessageEol;
3336            //realDualInfeasibilities = sumDualInfeasibilities_;
3337          } else if (numberChangedBounds<-1) {
3338            // bounds were flipped
3339            //gutsOfSolution(2);
3340          }
3341          if (stateDualColumn_>=0)
3342            stateDualColumn_++;
3343#ifndef TRY_ABC_GUS
3344          int returnCode=bounceTolerances(1);
3345          if (returnCode)
3346#endif // always go to primal
3347            problemStatus_=10;
3348          lastCleaned_ = numberIterations_;
3349#if CLP_CAN_HAVE_ZERO_OBJ
3350        } else {
3351          // no cost - skip checks
3352          problemStatus_=0;
3353        }
3354#endif
3355      } else {
3356        problemStatus_ = 0; // optimal
3357        if (lastCleaned_ < numberIterations_ && numberChanged_) {
3358          handler_->message(CLP_SIMPLEX_GIVINGUP, messages_)
3359            << CoinMessageEol;
3360        }
3361      }
3362    } else if (numberChangedBounds<-1) {
3363      // some flipped
3364      bounceTolerances(2);
3365#if ABC_NORMAL_DEBUG>0
3366      printf("numberChangedBounds %d numberAtFakeBound %d at line %d in file %s\n",
3367             numberChangedBounds,numberAtFakeBound(),__LINE__,__FILE__);
3368#endif
3369    } else if (numberAtFakeBound()) {
3370      handler_->message(CLP_DUAL_BOUNDS, messages_)
3371        << currentDualBound_
3372        << CoinMessageEol;
3373#if ABC_NORMAL_DEBUG>0
3374      printf("changing bounds\n");
3375      printf("numberChangedBounds %d numberAtFakeBound %d at line %d in file %s\n",
3376             numberChangedBounds,numberAtFakeBound(),__LINE__,__FILE__);
3377#endif
3378      if (!numberDualInfeasibilities_&&!numberPrimalInfeasibilities_) {
3379        // check unbounded
3380        // find a variable with bad dj
3381        int iChosen = -1;
3382        if ((specialOptions_ & 0x03000000) == 0) {
3383          double largest = 100.0 * primalTolerance_;
3384          for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
3385            double djValue = abcDj_[iSequence];
3386            double originalLo = lowerSaved_[iSequence];
3387            double originalUp = upperSaved_[iSequence];
3388            if (fabs(djValue) > fabs(largest)) {
3389              if (getInternalStatus(iSequence) != basic) {
3390                if (djValue > 0 && originalLo < -1.0e20) {
3391                  if (djValue > fabs(largest)) {
3392                    largest = djValue;
3393                    iChosen = iSequence;
3394                  }
3395                } else if (djValue < 0 && originalUp > 1.0e20) {
3396                  if (-djValue > fabs(largest)) {
3397                    largest = djValue;
3398                    iChosen = iSequence;
3399                  }
3400                }
3401              }
3402            }
3403          }
3404        }
3405        if (iChosen >= 0) {
3406          int iVector=getAvailableArray();
3407          unpack(usefulArray_[iVector],iChosen);
3408          //double changeCost;
3409          problemStatus_ = checkUnbounded(usefulArray_[iVector], 0.0/*changeCost*/);
3410          usefulArray_[iVector].clear();
3411          setAvailableArray(iVector);
3412        }
3413      }
3414      gutsOfSolution(3);
3415      //bounceTolerances(2);
3416      //realDualInfeasibilities = sumDualInfeasibilities_;
3417    }
3418    // unflag all variables (we may want to wait a bit?)
3419    if (unflagVariables) {
3420      int numberFlagged = numberFlagged_;
3421      numberFlagged_=0;
3422      for (int iRow = 0; iRow < numberRows_; iRow++) {
3423        int iPivot = abcPivotVariable_[iRow];
3424        if (flagged(iPivot)) {
3425          clearFlagged(iPivot);
3426        }
3427      }
3428#if ABC_NORMAL_DEBUG>3
3429      if (numberFlagged) {
3430        printf("unflagging %d variables - tentativeStatus %d probStat %d ninf %d nopt %d\n", numberFlagged, tentativeStatus,
3431               problemStatus_, numberPrimalInfeasibilities_,
3432               numberTimesOptimal_);
3433      }
3434#endif
3435      unflagVariables = numberFlagged > 0;
3436      if (numberFlagged && !numberPivots) {
3437        /* looks like trouble as we have not done any iterations.
3438           Try changing pivot tolerance then give it a few goes and give up */
3439        if (abcFactorization_->pivotTolerance() < 0.9) {
3440          abcFactorization_->pivotTolerance(0.99);
3441          problemStatus_ = -1;
3442        } else if (numberTimesOptimal_ < 3) {
3443          numberTimesOptimal_++;
3444          problemStatus_ = -1;
3445        } else {
3446          unflagVariables = false;
3447          //secondaryStatus_ = 1; // and say probably infeasible
3448          if ((moreSpecialOptions_ & 256) == 0||(abcState_&CLP_ABC_BEEN_FEASIBLE)!=0) {
3449            // try primal
3450            problemStatus_ = 10;
3451          } else {
3452            // almost certainly infeasible
3453            problemStatus_ = 1;
3454          }
3455#if ABC_NORMAL_DEBUG>3
3456          printf("returning at %d\n", __LINE__);
3457#endif
3458        }
3459      }
3460    }
3461  }
3462  if (problemStatus_ < 0) {
3463    saveGoodStatus();
3464    if (weightsSaved) {
3465      // restore weights (if saved) - also recompute infeasibility list
3466      abcDualRowPivot_->saveWeights(this, weightsSaved);
3467    } else {
3468      abcDualRowPivot_->recomputeInfeasibilities();
3469    }
3470  }
3471  // see if cutoff reached
3472  checkCutoff(false);
3473  // If we are in trouble and in branch and bound give up
3474  if ((specialOptions_ & 1024) != 0) {
3475    int looksBad = 0;
3476    if (largestPrimalError_ * largestDualError_ > 1.0e2) {
3477      looksBad = 1;
3478    } else if (largestPrimalError_ > 1.0e-2
3479               && rawObjectiveValue_ > CoinMin(1.0e15, 1.0e3 * limit)) {
3480      looksBad = 2;
3481    }
3482    if (looksBad) {
3483      if (abcFactorization_->pivotTolerance() < 0.9) {
3484        // up tolerance
3485        abcFactorization_->pivotTolerance(CoinMin(abcFactorization_->pivotTolerance() * 1.05 + 0.02, 0.91));
3486      } else if (numberIterations_ > 10000) {
3487#if ABC_NORMAL_DEBUG>0
3488        if (handler_->logLevel() > 2)
3489          printf("bad dual - saying infeasible %d\n", looksBad);
3490#endif
3491        problemStatus_ = 1;
3492        secondaryStatus_ = 1; // and say was on cutoff
3493      } else if (largestPrimalError_ > 1.0e5) {
3494#if ABC_NORMAL_DEBUG>0
3495        if (handler_->logLevel() > 2)
3496          printf("bad dual - going to primal %d %g\n", looksBad, largestPrimalError_);
3497#endif
3498        allSlackBasis();
3499        problemStatus_ = 10;
3500      }
3501    }
3502  }
3503  //if (problemStatus_ < 0 && !changeMade_) {
3504  //problemStatus_ = 4; // unknown
3505  //}
3506  lastGoodIteration_ = numberIterations_;
3507  if (numberIterations_ > lastBadIteration_ + 200) {
3508    moreSpecialOptions_ &= ~16; // clear check accuracy flag
3509    if (numberFlagged_) {
3510      numberFlagged_=0;
3511      for (int iRow = 0; iRow < numberRows_; iRow++) {
3512        int iPivot = abcPivotVariable_[iRow];
3513        if (flagged(iPivot)) {
3514          clearFlagged(iPivot);
3515        }
3516      }
3517    }
3518  }
3519  if (problemStatus_ < 0&&ignoreStuff!=2) {
3520    //sumDualInfeasibilities_ = realDualInfeasibilities; // back to say be careful
3521    if (sumDualInfeasibilities_) {
3522      numberDualInfeasibilities_ = 1;
3523    } else if (!numberPrimalInfeasibilities_) {
3524      // Looks good
3525      problemStatus_=0;
3526    }
3527  }
3528  double thisObj = abcProgress_.lastObjective(0);
3529  double lastObj = abcProgress_.lastObjective(1);
3530  if (lastObj > thisObj + 1.0e-4 * CoinMax(fabs(thisObj), fabs(lastObj)) + 1.0e-4&&
3531      lastFirstFree_<0) {
3532    //printf("BAD - objective backwards\n");
3533    //bounceTolerances(3);
3534    numberBlackMarks+=3;
3535  }
3536  if (numberBlackMarks>0&&numberIterations_>baseIteration_) {
3537    // be very careful
3538    int maxFactor = abcFactorization_->maximumPivots();
3539    if (maxFactor > 10) {
3540      if ((stateOfProblem_&PESSIMISTIC)==0||true) {
3541        if (largestDualError_>10.0*CoinMax(lastDualError_,1.0e-6)
3542            ||largestPrimalError_>10.0*CoinMax(lastPrimalError_,1.0e-6))
3543          maxFactor = CoinMin(maxFactor,20);
3544        forceFactorization_ = CoinMax(1, (maxFactor >> 1));
3545      } else if (numberBlackMarks>2) {
3546        forceFactorization_ = CoinMax(1, (forceFactorization_ >> 1));
3547      }
3548    }
3549    stateOfProblem_ |= PESSIMISTIC;
3550    if (numberBlackMarks>=10)
3551      forceFactorization_ = 1;
3552  }
3553  lastPrimalError_=largestPrimalError_;
3554  lastDualError_=largestDualError_;
3555  lastFirstFree_=firstFree_;
3556  if (ignoreStuff==2) {
3557    // in values pass
3558    lastFirstFree_=-1;
3559    // firstFree_=-1;
3560  }
3561  if (alphaAccuracy_ > 0.0)
3562    alphaAccuracy_ = 1.0;
3563  // If we are stopping - use plausible objective
3564  // Maybe only in fast dual
3565  if (problemStatus_ > 2)
3566    rawObjectiveValue_ = approximateObjective;
3567  if (problemStatus_ == 1 && (progressFlag_&8) != 0 &&
3568      fabs(rawObjectiveValue_) > 1.0e10 )
3569    problemStatus_ = 10; // infeasible - but has looked feasible
3570  checkMoveBack(true);
3571#if 0
3572  {
3573    double * array = usefulArray_[arrayForTableauRow_].denseVector();
3574    for (int iPivot = 0; iPivot < numberRows_; iPivot++) {
3575      int iSequence = abcPivotVariable_[iPivot];
3576      unpack(usefulArray_[arrayForTableauRow_],iSequence);
3577      abcFactorization_->updateColumn(usefulArray_[arrayForTableauRow_]);
3578      assert (fabs(array[iPivot] - 1.0) < 1.0e-4);
3579      array[iPivot] = 0.0;
3580      for (int i = 0; i < numberRows_; i++)
3581        assert (fabs(array[i]) < 1.0e-4);
3582      usefulArray_[arrayForTableauRow_].clear();
3583    }
3584  }
3585#endif
3586#if 0
3587  {
3588    // save status, cost and solution
3589    unsigned char * saveStatus = CoinCopyOfArray(internalStatus_,numberTotal_);
3590    double * saveSolution = CoinCopyOfArray(abcSolution_,numberTotal_);
3591    double * saveCost = CoinCopyOfArray(abcCost_,numberTotal_);
3592    printf("TestA sumDual %g sumPrimal %g\n",sumDualInfeasibilities_,sumPrimalInfeasibilities_);
3593    if (perturbation_>100)
3594      CoinAbcMemcpy(abcCost_,costSaved_,numberTotal_);
3595    else
3596      CoinAbcMemcpy(abcCost_,abcPerturbation_,numberTotal_);
3597    moveToBasic();
3598    gutsOfSolution(3);
3599    for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
3600      bool changed=fabs(saveCost[iSequence]-abcCost_[iSequence])>1.0e-10;
3601      bool bad=false;
3602      int iStatus=internalStatus_[iSequence]&7;
3603      if (getInternalStatus(iSequence) != basic && !flagged(iSequence)) {
3604        // not basic
3605        double distanceUp = abcUpper_[iSequence] - abcSolution_[iSequence];
3606        double distanceDown = abcSolution_[iSequence] - abcLower_[iSequence];
3607        double value = abcDj_[iSequence];
3608        if (distanceUp > primalTolerance_) {
3609          // Check if "free"
3610          if (distanceDown <= primalTolerance_) {
3611            // should not be negative
3612            if (value < 0.0) {
3613              value = - value;
3614              if (value > currentDualTolerance_) {
3615                bad=true;
3616              }
3617            }
3618          } else {
3619            // free
3620            value=fabs(value);
3621            if (value > currentDualTolerance_) {
3622              bad=true;
3623            }
3624          }
3625        } else if (distanceDown > primalTolerance_) {
3626          // should not be positive
3627          if (value > 0.0) {
3628            if (value > currentDualTolerance_) {
3629              bad=true;
3630            }
3631          }
3632        }
3633      }
3634      if (changed||bad)
3635        printf("seq %d status %d current cost %g original %g dj %g %s\n",
3636               iSequence,iStatus,saveCost[iSequence],abcCost_[iSequence],
3637               abcDj_[iSequence],bad ? "(bad)" : "");
3638    }
3639    printf("TestB sumDual %g sumPrimal %g\n",sumDualInfeasibilities_,sumPrimalInfeasibilities_);
3640    if (numberDualInfeasibilities_) {
3641      makeNonFreeVariablesDualFeasible();
3642      moveToBasic();
3643      gutsOfSolution(3);
3644      printf("TestC sumDual %g sumPrimal %g\n",sumDualInfeasibilities_,sumPrimalInfeasibilities_);
3645    }
3646    CoinAbcMemcpy(internalStatus_,saveStatus,numberTotal_);
3647    CoinAbcMemcpy(abcSolution_,saveSolution,numberTotal_);
3648    CoinAbcMemcpy(abcCost_,saveCost,numberTotal_);
3649    delete [] saveStatus;
3650    delete [] saveSolution;
3651    delete [] saveCost;
3652    moveToBasic();
3653    gutsOfSolution(3);
3654    printf("TestD sumDual %g sumPrimal %g\n",sumDualInfeasibilities_,sumPrimalInfeasibilities_);
3655  }
3656#endif
3657#if 0
3658  {
3659    const double multiplier[] = { 1.0, -1.0}; // sign?????????????
3660    double dualT = - currentDualTolerance_;
3661    double largestCostDifference=0.1*currentDualTolerance_;
3662    int numberCostDifference=0;
3663    double largestBadDj=0.0;
3664    int numberBadDj=0;
3665    double sumBadDj=0.0;
3666    double sumCanMoveBad=0.0;
3667    int numberCanMoveBad=0;
3668    double sumCanMoveOthers=0.0;
3669    int numberCanMoveOthers=0;
3670    for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
3671      double costDifference=abcPerturbation_[iSequence]-abcCost_[iSequence];
3672      if (fabs(costDifference)>0.1*currentDualTolerance_) {
3673        if (fabs(costDifference)>largestCostDifference)
3674          largestCostDifference=fabs(costDifference);
3675        numberCostDifference++;
3676      }
3677      unsigned char iStatus = internalStatus_[iSequence] & 7;
3678      if (iStatus<4) {
3679        double mult = multiplier[iStatus];
3680        double dj = abcDj_[iSequence] * mult;
3681        if (dj < dualT) {
3682          sumBadDj+=dualT-dj;
3683          numberBadDj++;
3684          if (dualT-dj>largestBadDj)
3685            largestBadDj=dualT-dj;
3686          if (costDifference*mult>0.0) {
3687            sumCanMoveBad+=costDifference*mult;
3688            numberCanMoveBad++;
3689            //dj+=costDifference*mult;
3690            abcDj_[iSequence]+=costDifference;
3691            abcCost_[iSequence]=abcPerturbation_[iSequence];
3692          }
3693        } else if (costDifference*mult>0.0) {
3694          sumCanMoveOthers+=costDifference*mult;
3695          numberCanMoveOthers++;
3696          //dj+=costDifference*mult;
3697          abcDj_[iSequence]+=costDifference;
3698          abcCost_[iSequence]=abcPerturbation_[iSequence];
3699        }
3700      }
3701    }
3702    if (numberCostDifference+numberBadDj+numberCanMoveBad+numberCanMoveOthers) {
3703      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,
3704             numberBadDj, sumBadDj,sumCanMoveBad=0.0,numberCanMoveBad,
3705             sumCanMoveOthers,numberCanMoveOthers,numberAtFakeBound());
3706    }
3707  }
3708#endif
3709}
3710// see if cutoff reached
3711bool
3712AbcSimplexDual::checkCutoff(bool computeObjective)
3713{
3714  double objValue=COIN_DBL_MAX;
3715  if (computeObjective)
3716    objValue = computeInternalObjectiveValue();
3717  else
3718    objValue = minimizationObjectiveValue();
3719  double limit = 0.0;
3720  getDblParam(ClpDualObjectiveLimit, limit);
3721  // Looks as if numberAtFakeBound() being computed too often
3722  if(fabs(limit) < 1.0e30 && objValue >
3723     limit && !numberAtFakeBound()) {
3724    bool looksInfeasible = !numberDualInfeasibilities_;
3725    if (objValue > limit + fabs(0.1 * limit) + 1.0e2 * sumDualInfeasibilities_ + 1.0e4 &&
3726        sumDualInfeasibilities_ < largestDualError_ && numberIterations_ > 0.5 * numberRows_ + 1000)
3727      looksInfeasible = true;
3728    if (looksInfeasible&&!computeObjective)
3729      objValue = computeInternalObjectiveValue();
3730    if (looksInfeasible) {
3731      // Even if not perturbed internal costs may have changed
3732      // be careful
3733      if(objValue > limit) {
3734        problemStatus_ = 1;
3735        secondaryStatus_ = 1; // and say was on cutoff
3736        return true;
3737      }
3738    }
3739  }
3740  return false;
3741}
3742#define CLEANUP_DJS 0
3743/* While updateDualsInDual sees what effect is of flip
3744   this does actual flipping.
3745   If change >0.0 then value in array >0.0 => from lower to upper
3746   returns 3 if skip this iteration and re-factorize
3747*/
3748int
3749AbcSimplexDual::flipBounds()
3750{
3751  CoinIndexedVector & array = usefulArray_[arrayForFlipBounds_];
3752  CoinIndexedVector & output = usefulArray_[arrayForFlipRhs_];
3753  output.clear();
3754  int numberFlipped=array.getNumElements();
3755  if (numberFlipped) {
3756#if CLP_CAN_HAVE_ZERO_OBJ>1
3757    if ((specialOptions_&2097152)==0) {
3758#endif
3759      // do actual flips
3760      // do in parallel once got spare space in factorization
3761#if ABC_PARALLEL==2
3762#endif
3763      int number=array.getNumElements();
3764      const int *  COIN_RESTRICT which=array.getIndices();
3765      double *  COIN_RESTRICT work = array.denseVector();
3766     
3767      double *  COIN_RESTRICT solution = abcSolution_;
3768      const double *  COIN_RESTRICT lower = abcLower_;
3769      const double *  COIN_RESTRICT upper = abcUpper_;
3770#if CLEANUP_DJS
3771      double *  COIN_RESTRICT cost = abcCost_;
3772      double *  COIN_RESTRICT dj = abcDj_;
3773      // see if we can clean up djs
3774      // may have perturbation
3775      const double * COIN_RESTRICT originalCost = abcPerturbation_;
3776#else
3777      const double *  COIN_RESTRICT cost = abcCost_;
3778      const double *  COIN_RESTRICT dj = abcDj_;
3779#endif
3780      const double multiplier[] = { 1.0, -1.0};
3781      for (int i = 0; i < number; i++) {
3782        double value=work[i]*theta_;
3783        work[i]=0.0;
3784        int iSequence = which[i];
3785        unsigned char iStatus = internalStatus_[iSequence]&3;
3786        assert ((internalStatus_[iSequence]&7)==0||
3787                (internalStatus_[iSequence]&7)==1);
3788        double mult = multiplier[iStatus];
3789        double djValue = dj[iSequence] * mult-value;
3790        //double perturbedTolerance = abcPerturbation_[iSequence];
3791        if (djValue<-currentDualTolerance_) {
3792          // might just happen - if so just skip
3793          assert (iSequence!=sequenceIn_);
3794          double movement=(upper[iSequence]-lower[iSequence])*mult;
3795          if (iStatus==0) {
3796            // at lower bound
3797            // to upper bound
3798            setInternalStatus(iSequence, atUpperBound);
3799            solution[iSequence] = upper[iSequence];
3800#if CLEANUP_DJS
3801            if (cost[iSequence]<originalCost[iSequence]) {
3802              double difference = CoinMax(cost[iSequence]-originalCost[iSequence],djValue);
3803              dj[iSequence]-=difference;
3804              cost[iSequence]-=difference;
3805            }
3806#endif
3807          } else {
3808            // at upper bound
3809            // to lower bound
3810            setInternalStatus(iSequence, atLowerBound);
3811            solution[iSequence] = lower[iSequence];
3812#if CLEANUP_DJS
3813            if (cost[iSequence]>originalCost[iSequence]) {
3814              double difference = CoinMin(cost[iSequence]-originalCost[iSequence],-djValue);
3815              dj[iSequence]-=difference;
3816              cost[iSequence]-=difference;
3817            }
3818#endif
3819          }
3820          abcMatrix_->add(output, iSequence, movement);
3821          objectiveChange_ += cost[iSequence]*movement;
3822        }
3823      }
3824      array.setNumElements(0);
3825      // double check something flipped
3826      if (output.getNumElements()) {
3827#if ABC_PARALLEL==0
3828        abcFactorization_->updateColumn(output);
3829#else
3830        assert (FACTOR_CPU>2);
3831        abcFactorization_->updateColumnCpu(output,1);
3832#endif
3833        abcDualRowPivot_->updatePrimalSolution(output,1.0);
3834      }
3835#if CLP_CAN_HAVE_ZERO_OBJ>1
3836    } else {
3837      clearArrays(3);
3838    }
3839#endif
3840  }
3841  // recompute dualOut_
3842  valueOut_ = solutionBasic_[pivotRow_];
3843  if (directionOut_ < 0) {
3844    dualOut_ = valueOut_ - upperOut_;
3845  } else {
3846    dualOut_ = lowerOut_ - valueOut_;
3847  }
3848#if 0
3849  // amount primal will move
3850  movement_ = -dualOut_ * directionOut_ / alpha_;
3851  double movementOld = oldDualOut * directionOut_ / alpha_;
3852  // so objective should increase by fabs(dj)*movement_
3853  // but we already have objective change - so check will be good
3854  if (objectiveChange_ + fabs(movementOld * dualIn_) < -CoinMax(1.0e-5, 1.0e-12 * fabs(rawObjectiveValue_))) {
3855#if ABC_NORMAL_DEBUG>3
3856    if (handler_->logLevel() & 32)
3857      printf("movement %g, movementOld %g swap change %g, rest %g  * %g\n",
3858             objectiveChange_ + fabs(movement_ * dualIn_),
3859             movementOld,objectiveChange_, movement_, dualIn_);
3860#endif
3861  }
3862  if (0) {
3863    if(abcFactorization_->pivots()) {
3864      // going backwards - factorize
3865      problemStatus_ = -2; // factorize now
3866      stateOfIteration_=2;
3867    }
3868  }
3869#endif
3870  return numberFlipped;
3871}
3872/* Undo a flip
3873 */
3874void 
3875AbcSimplexDual::flipBack(int number)
3876{
3877  CoinIndexedVector & array = usefulArray_[arrayForFlipBounds_];
3878  const int *  COIN_RESTRICT which=array.getIndices();
3879  double *  COIN_RESTRICT solution = abcSolution_;
3880  const double *  COIN_RESTRICT lower = abcLower_;
3881  const double *  COIN_RESTRICT upper = abcUpper_;
3882  for (int i = 0; i < number; i++) {
3883    int iSequence = which[i];
3884    unsigned char iStatus = internalStatus_[iSequence]&3;
3885    assert ((internalStatus_[iSequence]&7)==0||
3886            (internalStatus_[iSequence]&7)==1);
3887    if (iStatus==0) {
3888      // at lower bound
3889      // to upper bound
3890      setInternalStatus(iSequence, atUpperBound);
3891      solution[iSequence] = upper[iSequence];
3892    } else {
3893      // at upper bound
3894      // to lower bound
3895      setInternalStatus(iSequence, atLowerBound);
3896      solution[iSequence] = lower[iSequence];
3897    }
3898  }
3899}
3900// Restores bound to original bound
3901void
3902AbcSimplexDual::originalBound( int iSequence)
3903{
3904  if (getFakeBound(iSequence) != noFake) {
3905    numberFake_--;;
3906    setFakeBound(iSequence, noFake);
3907    abcLower_[iSequence] = lowerSaved_[iSequence];
3908    abcUpper_[iSequence] = upperSaved_[iSequence];
3909  }
3910}
3911/* As changeBounds but just changes new bounds for a single variable.
3912   Returns true if change */
3913bool
3914AbcSimplexDual::changeBound( int iSequence)
3915{
3916  // old values
3917  double oldLower = abcLower_[iSequence];
3918  double oldUpper = abcUpper_[iSequence];
3919  double value = abcSolution_[iSequence];
3920  bool modified = false;
3921  // original values
3922  double lowerValue = lowerSaved_[iSequence];
3923  double upperValue = upperSaved_[iSequence];
3924  setFakeBound(iSequence,noFake);
3925  if (value == oldLower) {
3926    if (upperValue > oldLower + currentDualBound_) {
3927      abcUpper_[iSequence] = oldLower + currentDualBound_;
3928      setFakeBound(iSequence, upperFake);
3929      assert (abcLower_[iSequence]==lowerSaved_[iSequence]);
3930      assert (abcUpper_[iSequence]<upperSaved_[iSequence]);
3931      modified = true;
3932      numberFake_++;
3933    }
3934  } else if (value == oldUpper) {
3935    if (lowerValue < oldUpper - currentDualBound_) {
3936      abcLower_[iSequence] = oldUpper - currentDualBound_;
3937      setFakeBound(iSequence, lowerFake);
3938      assert (abcLower_[iSequence]>lowerSaved_[iSequence]);
3939      assert (abcUpper_[iSequence]==upperSaved_[iSequence]);
3940      modified = true;
3941      numberFake_++;
3942    }
3943  } else {
3944    assert(value == oldLower || value == oldUpper);
3945  }
3946  return modified;
3947}
3948/* Checks number of variables at fake bounds.  This is used by fastDual
3949   so can exit gracefully before end */
3950int
3951AbcSimplexDual::numberAtFakeBound()
3952{
3953  int numberFake = 0;
3954
3955  if (numberFake_) {
3956    for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
3957      FakeBound bound = getFakeBound(iSequence);
3958      switch(getInternalStatus(iSequence)) {
3959       
3960      case basic:
3961        break;
3962      case isFree:
3963      case superBasic:
3964      case AbcSimplex::isFixed:
3965        //setFakeBound (iSequence, noFake);
3966        break;
3967      case atUpperBound:
3968        if (bound == upperFake || bound == bothFake)
3969          numberFake++;
3970        break;
3971      case atLowerBound:
3972        if (bound == lowerFake || bound == bothFake)
3973          numberFake++;
3974        break;
3975      }
3976    }
3977  }
3978  return numberFake;
3979}
3980int
3981AbcSimplexDual::resetFakeBounds(int type)
3982{
3983  if (type == 0||type==1) {
3984    // put back original bounds and then check
3985    //CoinAbcMemcpy(abcLower_,lowerSaved_,numberTotal_);
3986    //CoinAbcMemcpy(abcUpper_,upperSaved_,numberTotal_);
3987    double dummyChangeCost = 0.0;
3988    return changeBounds(3+type, dummyChangeCost);
3989  } else if (type < 0) {
3990#ifndef NDEBUG
3991    // just check
3992    int iSequence;
3993    int nFake = 0;
3994    int nErrors = 0;
3995    int nSuperBasic = 0;
3996    int nWarnings = 0;
3997    for (iSequence = 0; iSequence < numberTotal_; iSequence++) {
3998      FakeBound fakeStatus = getFakeBound(iSequence);
3999      Status status = getInternalStatus(iSequence);
4000      bool isFake = false;
4001#ifdef CLP_INVESTIGATE
4002      char RC = 'C';
4003      int jSequence = iSequence;
4004      if (jSequence >= numberColumns_) {
4005        RC = 'R';
4006        jSequence -= numberColumns_;
4007      }
4008#endif
4009      double lowerValue = lowerSaved_[iSequence];
4010      double upperValue = upperSaved_[iSequence];
4011      double value = abcSolution_[iSequence];
4012      CoinRelFltEq equal;
4013      if (status == atUpperBound ||
4014          status == atLowerBound) {
4015        if (fakeStatus == AbcSimplexDual::upperFake) {
4016          if(!equal(abcUpper_[iSequence], (lowerValue + currentDualBound_)) ||
4017             !(equal(abcUpper_[iSequence], value) ||
4018               equal(abcLower_[iSequence], value))) {
4019            nErrors++;
4020#ifdef CLP_INVESTIGATE
4021            printf("** upperFake %c%d %g <= %g <= %g true %g, %g\n",
4022                   RC, jSequence, abcLower_[iSequence], abcSolution_[iSequence],
4023                   abcUpper_[iSequence], lowerValue, upperValue);
4024#endif
4025          }
4026          isFake = true;;
4027        } else if (fakeStatus == AbcSimplexDual::lowerFake) {
4028          if(!equal(abcLower_[iSequence], (upperValue - currentDualBound_)) ||
4029             !(equal(abcUpper_[iSequence], value) ||
4030               equal(abcLower_[iSequence], value))) {
4031            nErrors++;
4032#ifdef CLP_INVESTIGATE
4033            printf("** lowerFake %c%d %g <= %g <= %g true %g, %g\n",
4034                   RC, jSequence, abcLower_[iSequence], abcSolution_[iSequence],
4035                   abcUpper_[iSequence], lowerValue, upperValue);
4036#endif
4037          }
4038          isFake = true;;
4039        } else if (fakeStatus == AbcSimplexDual::bothFake) {
4040          nWarnings++;
4041#ifdef CLP_INVESTIGATE
4042          printf("** %d at bothFake?\n", iSequence);
4043#endif
4044        } else if (abcUpper_[iSequence] - abcLower_[iSequence] > 2.0 * currentDualBound_) {
4045          nErrors++;
4046#ifdef CLP_INVESTIGATE
4047          printf("** noFake! %c%d %g <= %g <= %g true %g, %g\n",
4048                 RC, jSequence, abcLower_[iSequence], abcSolution_[iSequence],
4049                 abcUpper_[iSequence], lowerValue, upperValue);
4050#endif
4051        }
4052      } else if (status == superBasic || status == isFree) {
4053        nSuperBasic++;
4054        //printf("** free or superbasic %c%d %g <= %g <= %g true %g, %g - status %d\n",
4055        //     RC,jSequence,abcLower_[iSequence],abcSolution_[iSequence],
4056        //     abcUpper_[iSequence],lowerValue,upperValue,status);
4057      } else if (status == basic) {
4058        bool odd = false;
4059        if (!equal(abcLower_[iSequence], lowerValue))
4060          odd = true;
4061        if (!equal(abcUpper_[iSequence], upperValue))
4062          odd = true;
4063        if (odd) {
4064#ifdef CLP_INVESTIGATE
4065          printf("** basic %c%d %g <= %g <= %g true %g, %g\n",
4066                 RC, jSequence, abcLower_[iSequence], abcSolution_[iSequence],
4067                 abcUpper_[iSequence], lowerValue, upperValue);
4068#endif
4069          nWarnings++;
4070        }
4071      } else if (status == isFixed) {
4072        if (!equal(abcUpper_[iSequence], abcLower_[iSequence])) {
4073          nErrors++;
4074#ifdef CLP_INVESTIGATE
4075          printf("** fixed! %c%d %g <= %g <= %g true %g, %g\n",
4076                 RC, jSequence, abcLower_[iSequence], abcSolution_[iSequence],
4077                 abcUpper_[iSequence], lowerValue, upperValue);
4078#endif
4079        }
4080      }
4081      if (isFake) {
4082        nFake++;
4083      } else {
4084        if (fakeStatus != AbcSimplexDual::noFake) {
4085          nErrors++;
4086#ifdef CLP_INVESTIGATE
4087          printf("** bad fake status %c%d %d\n",
4088                 RC, jSequence, fakeStatus);
4089#endif
4090        }
4091      }
4092    }
4093    if (nFake != numberFake_) {
4094#ifdef CLP_INVESTIGATE
4095      printf("nfake %d numberFake %d\n", nFake, numberFake_);
4096#endif
4097      nErrors++;
4098    }
4099    if (nErrors || type <= -1000) {
4100#ifdef CLP_INVESTIGATE
4101      printf("%d errors, %d warnings, %d free/superbasic, %d fake\n",
4102             nErrors, nWarnings, nSuperBasic, numberFake_);
4103      printf("currentDualBound %g\n",
4104             currentDualBound_);
4105#endif
4106      if (type <= -1000) {
4107        iSequence = -type;
4108        iSequence -= 1000;
4109#ifdef CLP_INVESTIGATE
4110        char RC = 'C';
4111        int jSequence = iSequence;
4112        if (jSequence >= numberColumns_) {
4113          RC = 'R';
4114          jSequence -= numberColumns_;
4115        }
4116        double lowerValue = lowerSaved_[iSequence];
4117        double upperValue = upperSaved_[iSequence];
4118        printf("*** movement>1.0e30 for  %c%d %g <= %g <= %g true %g, %g - status %d\n",
4119               RC, jSequence, abcLower_[iSequence], abcSolution_[iSequence],
4120               abcUpper_[iSequence], lowerValue, upperValue, internalStatus_[iSequence]);
4121#endif
4122        assert (nErrors); // should have been picked up
4123      }
4124      assert (!nErrors);
4125    }
4126#endif
4127  } else {
4128    CoinAbcMemcpy(abcLower_,lowerSaved_,numberTotal_);
4129    CoinAbcMemcpy(abcUpper_,upperSaved_,numberTotal_);
4130    // reset using status
4131    numberFake_ = 0;
4132    for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
4133      FakeBound fakeStatus = getFakeBound(iSequence);
4134      if (fakeStatus != AbcSimplexDual::noFake) {
4135        Status status = getInternalStatus(iSequence);
4136        if (status == basic) {
4137          setFakeBound(iSequence, AbcSimplexDual::noFake);
4138          continue;
4139        }
4140        double lowerValue = abcLower_[iSequence];
4141        double upperValue = abcUpper_[iSequence];
4142        double value = abcSolution_[iSequence];
4143        numberFake_++;
4144        if (fakeStatus == AbcSimplexDual::upperFake) {
4145          abcUpper_[iSequence] = lowerValue + currentDualBound_;
4146          if (status == AbcSimplex::atLowerBound) {
4147            abcSolution_[iSequence] = lowerValue;
4148          } else if (status == AbcSimplex::atUpperBound) {
4149            abcSolution_[iSequence] = abcUpper_[iSequence];
4150          } else {
4151            abort();
4152          }
4153        } else if (fakeStatus == AbcSimplexDual::lowerFake) {
4154          abcLower_[iSequence] = upperValue - currentDualBound_;
4155          if (status == AbcSimplex::atLowerBound) {
4156            abcSolution_[iSequence] = abcLower_[iSequence];
4157          } else if (status == AbcSimplex::atUpperBound) {
4158            abcSolution_[iSequence] = upperValue;
4159          } else {
4160            abort();
4161          }
4162        } else {
4163          assert (fakeStatus == AbcSimplexDual::bothFake);
4164          if (status == AbcSimplex::atLowerBound) {
4165            abcLower_[iSequence] = value;
4166            abcUpper_[iSequence] = value + currentDualBound_;
4167          } else if (status == AbcSimplex::atUpperBound) {
4168            abcUpper_[iSequence] = value;
4169            abcLower_[iSequence] = value - currentDualBound_;
4170          } else if (status == AbcSimplex::isFree ||
4171                     status == AbcSimplex::superBasic) {
4172            abcLower_[iSequence] = value - 0.5 * currentDualBound_;
4173            abcUpper_[iSequence] = value + 0.5 * currentDualBound_;
4174          } else {
4175            abort();
4176          }
4177        }
4178      }
4179    }
4180  }
4181  return -1;
4182}
4183// maybe could JUST check against perturbationArray
4184void
4185AbcSimplex::checkDualSolutionPlusFake()
4186{
4187 
4188  sumDualInfeasibilities_ = 0.0;
4189  numberDualInfeasibilities_ = 0;
4190  sumFakeInfeasibilities_=0.0;
4191  int firstFreePrimal = -1;
4192  int firstFreeDual = -1;
4193  int numberSuperBasicWithDj = 0;
4194  bestPossibleImprovement_ = 0.0;
4195  // we can't really trust infeasibilities if there is dual error
4196  double error = CoinMin(1.0e-2, largestDualError_);
4197  // allow tolerance at least slightly bigger than standard
4198  double relaxedTolerance = currentDualTolerance_ +  error;
4199  // allow bigger tolerance for possible improvement
4200  double possTolerance = 5.0 * relaxedTolerance;
4201  sumOfRelaxedDualInfeasibilities_ = 0.0;
4202  //rawObjectiveValue_ -=sumNonBasicCosts_;
4203  sumNonBasicCosts_ = 0.0;
4204  const double *  COIN_RESTRICT abcLower = abcLower_;
4205  const double *  COIN_RESTRICT abcUpper = abcUpper_;
4206  const double *  COIN_RESTRICT abcSolution = abcSolution_;
4207#ifdef CLEANUP_DJS
4208  double *  COIN_RESTRICT abcDj = abcDj_;
4209  double *  COIN_RESTRICT abcCost = abcCost_;
4210  // see if we can clean up djs
4211  // may have perturbation
4212  const double * COIN_RESTRICT originalCost = abcPerturbation_;
4213#else
4214  const double *  COIN_RESTRICT abcDj = abcDj_;
4215  const double *  COIN_RESTRICT abcCost = abcCost_;
4216#endif
4217  //const double *  COIN_RESTRICT abcPerturbation = abcPerturbation_;
4218  int badDjSequence=-1;
4219  for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
4220    if (getInternalStatus(iSequence) != basic) {
4221      sumNonBasicCosts_ += abcSolution[iSequence]*abcCost[iSequence];
4222      // not basic
4223      double distanceUp = abcUpper[iSequence] - abcSolution[iSequence];
4224      double distanceDown = abcSolution[iSequence] - abcLower[iSequence];
4225      double value = abcDj[iSequence];
4226      if (distanceUp > primalTolerance_) {
4227        // Check if "free"
4228        if (distanceDown <= primalTolerance_) {
4229          // should not be negative
4230          if (value < 0.0) {
4231            value = - value;
4232            if (value > currentDualTolerance_) {
4233              bool treatAsSuperBasic=false;
4234              if (getInternalStatus(iSequence)==superBasic) {
4235#ifdef PAN
4236                if (fakeSuperBasic(iSequence)>=0) {
4237#endif
4238                  treatAsSuperBasic=true;
4239                  numberSuperBasicWithDj++;
4240                  if (firstFreeDual < 0)
4241                    firstFreeDual = iSequence;
4242#ifdef PAN
4243                }
4244#endif
4245              }
4246              if (!treatAsSuperBasic) {
4247                sumDualInfeasibilities_ += value - currentDualTolerance_;
4248                badDjSequence=iSequence;
4249                if (value > possTolerance)
4250                  bestPossibleImprovement_ += CoinMin(distanceUp, 1.0e10) * value;
4251                if (value > relaxedTolerance)
4252                  sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
4253                numberDualInfeasibilities_ ++;
4254              }
4255            }
4256#ifdef CLEANUP_DJS
4257          } else {
4258            // adjust cost etc
4259            if (abcCost[iSequence]>originalCost[iSequence]) {
4260              double difference = CoinMin(abcCost[iSequence]-originalCost[iSequence],value);
4261              //printf("Lb %d changing cost from %g to %g, dj from %g to %g\n",
4262              //     iSequence,abcCost[iSequence],abcCost[iSequence]-difference,
4263              //     abcDj[iSequence],abcDj[iSequence]-difference);
4264              abcDj[iSequence]=value-difference;
4265              abcCost[iSequence]-=difference;
4266            }
4267#endif
4268          }
4269        } else {
4270          // free
4271          value=fabs(value);
4272          if (value > 1.0 * relaxedTolerance) {
4273            numberSuperBasicWithDj++;
4274            if (firstFreeDual < 0)
4275              firstFreeDual = iSequence;
4276          }
4277          if (value > currentDualTolerance_) {
4278            sumDualInfeasibilities_ += value - currentDualTolerance_;
4279            if (value > possTolerance)
4280              bestPossibleImprovement_ += CoinMin(distanceUp, 1.0e10) * value;
4281            if (value > relaxedTolerance)
4282              sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
4283            numberDualInfeasibilities_ ++;
4284          }
4285          if (firstFreePrimal < 0)
4286            firstFreePrimal = iSequence;
4287        }
4288      } else if (distanceDown > primalTolerance_) {
4289        // should not be positive
4290        if (value > 0.0) {
4291          if (value > currentDualTolerance_) {
4292            bool treatAsSuperBasic=false;
4293            if (getInternalStatus(iSequence)==superBasic) {
4294#ifdef PAN
4295              if (fakeSuperBasic(iSequence)>=0) {
4296#endif
4297                treatAsSuperBasic=true;
4298                numberSuperBasicWithDj++;
4299                if (firstFreeDual < 0)
4300                  firstFreeDual = iSequence;
4301#ifdef PAN
4302              }
4303#endif
4304            }
4305            if (!treatAsSuperBasic) {
4306              sumDualInfeasibilities_ += value - currentDualTolerance_;
4307              badDjSequence=iSequence;
4308              if (value > possTolerance)
4309                bestPossibleImprovement_ += value * CoinMin(distanceDown, 1.0e10);
4310              if (value > relaxedTolerance)
4311                sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
4312              numberDualInfeasibilities_ ++;
4313            }
4314          }
4315#ifdef CLEANUP_DJS
4316        } else {
4317          // adjust cost etc
4318          if (abcCost[iSequence]<originalCost[iSequence]) {
4319            double difference = CoinMax(abcCost[iSequence]-originalCost[iSequence],value);
4320            //printf("Ub %d changing cost from %g to %g, dj from %g to %g\n",
4321            //       iSequence,abcCost[iSequence],abcCost[iSequence]-difference,
4322            //       abcDj[iSequence],abcDj[iSequence]-difference);
4323            abcDj[iSequence]=value-difference;
4324            abcCost[iSequence]-=difference;
4325          }
4326#endif
4327        }
4328      }
4329    }
4330  }
4331#ifdef CLEANUP_DJS
4332  if (algorithm_<0&&numberDualInfeasibilities_==1&&firstFreeDual<0&&sumDualInfeasibilities_<1.0e-1) {
4333    // clean bad one
4334    assert (badDjSequence>=0);
4335    abcCost[badDjSequence]-=abcDj_[badDjSequence];
4336    abcDj_[badDjSequence]=0.0;
4337  }
4338#endif
4339  numberDualInfeasibilitiesWithoutFree_ = numberDualInfeasibilities_-numberSuperBasicWithDj;
4340  numberFreeNonBasic_=numberSuperBasicWithDj;
4341  firstFree_=-1;
4342  if (algorithm_ < 0 && firstFreeDual >= 0) {
4343    // dual
4344    firstFree_ = firstFreeDual;
4345  } else if ((numberSuperBasicWithDj ||
4346              (abcProgress_.lastIterationNumber(0) <= 0))&&algorithm_>0) {
4347    firstFree_ = firstFreePrimal;
4348  }
4349}
4350// Perturbs problem
4351void 
4352AbcSimplexDual::perturb(double factor)
4353{
4354  const double * COIN_RESTRICT lower = lowerSaved_;
4355  const double * COIN_RESTRICT upper = upperSaved_;
4356  const double * COIN_RESTRICT perturbation = abcPerturbation_;
4357  double * COIN_RESTRICT cost = abcCost_;
4358  for (int i=0;i<numberTotal_;i++) {
4359    if (getInternalStatus(i)!=basic) {
4360      if (fabs(lower[i])<fabs(upper[i])) {
4361        cost[i]+=factor*perturbation[i];
4362      } else if (upper[i]!=COIN_DBL_MAX&&upper[i]!=lower[i]) {
4363        cost[i]-=factor*perturbation[i];
4364      } else {
4365        // free
4366        //cost[i]=costSaved[i];
4367      }
4368    }
4369  }
4370#ifdef CLEANUP_DJS
4371  // copy cost to perturbation
4372  CoinAbcMemcpy(abcPerturbation_,abcCost_,numberTotal_);
4373#endif
4374}
4375#if ABC_NORMAL_DEBUG>0
4376#define PERT_STATISTICS
4377#endif
4378#ifdef PERT_STATISTICS
4379static void breakdown(const char * name, int numberLook, const double * region)
4380{
4381     double range[] = {
4382          -COIN_DBL_MAX,
4383          -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,
4384          -1.0, -3.0e-1,
4385          -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,
4386          0.0,
4387          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,
4388          1.0, 3.0,
4389          1.0e1, 3.0e1, 1.0e2, 3.0e2, 1.0e3, 3.0e3, 1.0e4, 3.0e4, 1.0e5, 1.0e8, 1.0e11, 1.0e15,
4390          COIN_DBL_MAX
4391     };
4392     int nRanges = static_cast<int> (sizeof(range) / sizeof(double));
4393     int * number = new int[nRanges];
4394     memset(number, 0, nRanges * sizeof(int));
4395     int * numberExact = new int[nRanges];
4396     memset(numberExact, 0, nRanges * sizeof(int));
4397     int i;
4398     for ( i = 0; i < numberLook; i++) {
4399          double value = region[i];
4400          for (int j = 0; j < nRanges; j++) {
4401               if (value == range[j]) {
4402                    numberExact[j]++;
4403                    break;
4404               } else if (value < range[j]) {
4405                    number[j]++;
4406                    break;
4407               }
4408          }
4409     }
4410     printf("\n%s has %d entries\n", name, numberLook);
4411     for (i = 0; i < nRanges; i++) {
4412          if (number[i])
4413               printf("%d between %g and %g", number[i], range[i-1], range[i]);
4414          if (numberExact[i]) {
4415               if (number[i])
4416                    printf(", ");
4417               printf("%d exactly at %g", numberExact[i], range[i]);
4418          }
4419          if (number[i] + numberExact[i])
4420               printf("\n");
4421     }
4422     delete [] number;
4423     delete [] numberExact;
4424}
4425#endif
4426// Perturbs problem
4427void 
4428AbcSimplexDual::perturbB(double /*factorIn*/,int /*type*/)
4429{
4430  double overallMultiplier= (perturbation_==50||perturbation_>54) ? 1.0 : 0.1;
4431  // dual perturbation
4432  double perturbation = 1.0e-20;
4433  // maximum fraction of cost to perturb
4434  double maximumFraction = 1.0e-5;
4435  int maxLength = 0;
4436  int minLength = numberRows_;
4437  double averageCost = 0.0;
4438  int numberNonZero = 0;
4439  // See if we need to perturb
4440  double * COIN_RESTRICT sort = new double[numberColumns_];
4441  // Use objective BEFORE scaling ??
4442  const double * COIN_RESTRICT obj = objective();
4443  for (int i = 0; i < numberColumns_; i++) {
4444    double value = fabs(obj[i]);
4445    sort[i] = value;
4446    averageCost += fabs(abcCost_[i+maximumAbcNumberRows_]);
4447    if (value)
4448      numberNonZero++;
4449  }
4450  if (numberNonZero)
4451    averageCost /= static_cast<double> (numberNonZero);
4452  else
4453    averageCost = 1.0;
4454  // allow for cost scaling
4455  averageCost *= objectiveScale_;
4456  std::sort(sort, sort + numberColumns_);
4457  int number = 1;
4458  double last = sort[0];
4459  for (int i = 1; i < numberColumns_; i++) {
4460    if (last != sort[i])
4461      number++;
4462    last = sort[i];
4463  }
4464  delete [] sort;
4465#ifdef PERT_STATISTICS
4466  printf("%d different values\n",number);
4467#endif
4468#if 0
4469  printf("nnz %d percent %d", number, (number * 100) / numberColumns_);
4470  if (number * 4 > numberColumns_)
4471    printf(" - Would not perturb\n");
4472  else
4473    printf(" - Would perturb\n");
4474#endif
4475  // If small costs then more dangerous (better to scale)
4476  double numberMultiplier=1.25;
4477  if (averageCost<1.0e-3)
4478    numberMultiplier=20.0;
4479  else if (averageCost<1.0e-2)
4480    numberMultiplier=10.0;
4481  else if (averageCost<1.0e-1)
4482    numberMultiplier=3.0;
4483  else if (averageCost<1.0)
4484    numberMultiplier=2.0;
4485  numberMultiplier *= 2.0; // try
4486  if (number * numberMultiplier > numberColumns_||perturbation_>=100) {
4487    perturbation_=101;
4488#if ABC_NORMAL_DEBUG>0
4489    printf("Not perturbing - %d different costs, average %g\n",
4490           number,averageCost);
4491#endif
4492#ifdef CLEANUP_DJS
4493    // copy cost to perturbation
4494    CoinAbcMemcpy(abcPerturbation_,abcCost_,numberTotal_);
4495#endif
4496    return ; // good enough
4497  }
4498  const int * COIN_RESTRICT columnLength = matrix_->getVectorLengths()-maximumAbcNumberRows_;
4499  const double * COIN_RESTRICT lower = lowerSaved_;
4500  const double * COIN_RESTRICT upper = upperSaved_;
4501  for (int iSequence = maximumAbcNumberRows_; iSequence < numberTotal_; iSequence++) {
4502    if (lower[iSequence] < upper[iSequence]) {
4503      int length = columnLength[iSequence];
4504      if (length > 2) {
4505        maxLength = CoinMax(maxLength, length);
4506        minLength = CoinMin(minLength, length);
4507      }
4508    }
4509  }
4510  bool uniformChange = false;
4511  bool inCbcOrOther = (specialOptions_ & 0x03000000) != 0;
4512  double constantPerturbation = 10.0 * dualTolerance_;
4513#ifdef PERT_STATISTICS
4514  breakdown("Objective before perturbation", numberColumns_, abcCost_+numberRows_);
4515#endif
4516  //#define PERTURB_OLD_WAY
4517#ifndef PERTURB_OLD_WAY
4518#if 1
4519#ifdef HEAVY_PERTURBATION
4520    if (perturbation_==50)
4521      perturbation_=HEAVY_PERTURBATION;
4522#endif
4523  if (perturbation_ >= 50) {
4524    // Experiment
4525    // maximumFraction could be 1.0e-10 to 1.0
4526    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};
4527    int whichOne = CoinMin(perturbation_ - 51,10);
4528    if (whichOne<0) {
4529      if (numberRows_<5000)
4530        whichOne=2; // treat 50 as if 53
4531      else
4532        whichOne=2; // treat 50 as if 52
4533    }
4534    if (inCbcOrOther&&whichOne>0)
4535      whichOne--;
4536    maximumFraction = m[whichOne];
4537    if (whichOne>7)
4538      constantPerturbation *= 10.0; 
4539    //if (whichOne<2)
4540    //constantPerturbation = maximumFraction*1.0e-10 * dualTolerance_;
4541  } else if (inCbcOrOther) {
4542    maximumFraction = 1.0e-6;
4543  }
4544#endif
4545  double smallestNonZero = 1.0e100;
4546  bool modifyRowCosts=false;
4547  numberNonZero = 0;
4548  //perturbation = 1.0e-10;
4549  perturbation = CoinMax(1.0e-10,maximumFraction);
4550  double * COIN_RESTRICT cost = abcCost_;
4551  bool allSame = true;
4552  double lastValue = 0.0;
4553 
4554  for (int iRow = 0; iRow < numberRows_; iRow++) {
4555    double lo = lower[iRow];
4556    double up = upper[iRow];
4557    if (lo < up) {
4558      double value = fabs(cost[iRow]);
4559      perturbation = CoinMax(perturbation, value);
4560      if (value) {
4561        modifyRowCosts = true;
4562        smallestNonZero = CoinMin(smallestNonZero, value);
4563      }
4564    }
4565    if (lo && lo > -1.0e10) {
4566      numberNonZero++;
4567      lo = fabs(lo);
4568      if (!lastValue)
4569        lastValue = lo;
4570      else if (fabs(lo - lastValue) > 1.0e-7)
4571        allSame = false;
4572    }
4573    if (up && up < 1.0e10) {
4574      numberNonZero++;
4575      up = fabs(up);
4576      if (!lastValue)
4577        lastValue = up;
4578      else if (fabs(up - lastValue) > 1.0e-7)
4579        allSame = false;
4580    }
4581  }
4582  double lastValue2 = 0.0;
4583  for (int iSequence = maximumAbcNumberRows_; iSequence < numberTotal_; iSequence++) {
4584    double lo = lower[iSequence];
4585    double up = upper[iSequence];
4586    if (lo < up) {
4587      double value =
4588        fabs(cost[iSequence]);
4589      perturbation = CoinMax(perturbation, value);
4590      if (value) {
4591        smallestNonZero = CoinMin(smallestNonZero, value);
4592      }
4593    }
4594    if (lo && lo > -1.0e10) {
4595      //numberNonZero++;
4596      lo = fabs(lo);
4597      if (!lastValue2)
4598        lastValue2 = lo;
4599      else if (fabs(lo - lastValue2) > 1.0e-7)
4600        allSame = false;
4601    }
4602    if (up && up < 1.0e10) {
4603      //numberNonZero++;
4604      up = fabs(up);
4605      if (!lastValue2)
4606        lastValue2 = up;
4607      else if (fabs(up - lastValue2) > 1.0e-7)
4608        allSame = false;
4609    }
4610  }
4611  if (allSame) {
4612    // Check elements
4613    double smallestNegative;
4614    double largestNegative;
4615    double smallestPositive;
4616    double largestPositive;
4617    matrix_->rangeOfElements(smallestNegative, largestNegative,
4618                             smallestPositive, largestPositive);
4619    if (smallestNegative == largestNegative &&
4620        smallestPositive == largestPositive) {
4621      // Really hit perturbation
4622      double adjust = CoinMin(100.0 * maximumFraction, 1.0e-3 * CoinMax(lastValue, lastValue2));
4623      maximumFraction = CoinMax(adjust, maximumFraction);
4624    }
4625  }
4626  perturbation = CoinMin(perturbation, smallestNonZero / maximumFraction);
4627  double largestZero = 0.0;
4628  double largest = 0.0;
4629  double largestPerCent = 0.0;
4630  // modify costs
4631  bool printOut = (handler_->logLevel() == 63);
4632  printOut = false;
4633  //assert (!modifyRowCosts);
4634  modifyRowCosts = false;
4635  const double * COIN_RESTRICT perturbationArray = perturbationSaved_;
4636  if (modifyRowCosts) {
4637    for (int iRow = 0; iRow < numberRows_; iRow++) {
4638      if (lower[iRow] < upper[iRow]) {
4639        double value = perturbation;
4640        double currentValue = cost[iRow];
4641        value = CoinMin(value, maximumFraction * (fabs(currentValue) + 1.0e-1 * perturbation + 1.0e-3));
4642        double perturbationValue=2.0*(perturbationArray[iRow]-0.5); // for now back to normal random
4643        if (lower[iRow] > -largeValue_) {
4644          if (fabs(lower[iRow]) < fabs(upper[iRow]))
4645            value *= perturbationValue;
4646          else
4647            value *= -perturbationValue;
4648        } else if (upper[iRow] < largeValue_) {
4649          value *= -perturbationValue;
4650        } else {
4651          value = 0.0;
4652        }
4653        if (currentValue) {
4654          largest = CoinMax(largest, fabs(value));
4655          if (fabs(value) > fabs(currentValue)*largestPerCent)
4656            largestPerCent = fabs(value / currentValue);
4657        } else {
4658          largestZero = CoinMax(largestZero, fabs(value));
4659        }
4660        if (printOut)
4661          printf("row %d cost %g change %g\n", iRow, cost[iRow], value);
4662        cost[iRow] += value;
4663      }
4664    }
4665  }
4666
4667  //double extraWeight = 10.0;
4668  // 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};
4669  // 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};
4670  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};
4671  //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};
4672  // Scale back if wanted
4673  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};
4674  if (constantPerturbation < 99.0 * dualTolerance_&&false) {
4675    perturbation *= 0.1;
4676    //extraWeight = 0.5;
4677    memcpy(weight, weight2, sizeof(weight2));
4678  }
4679  // adjust weights if all columns long
4680  double factor = 1.0;
4681  if (maxLength) {
4682    factor = 3.0 / static_cast<double> (minLength);
4683  }
4684  // Make variables with more elements more expensive
4685  const double m1 = 0.5;
4686  double smallestAllowed = overallMultiplier*CoinMin(1.0e-2 * dualTolerance_, maximumFraction);
4687  double largestAllowed = overallMultiplier*CoinMax(1.0e3 * dualTolerance_, maximumFraction * averageCost);
4688  // smaller if in BAB
4689  //if (inCbcOrOther)
4690  //largestAllowed=CoinMin(largestAllowed,1.0e-5);
4691  //smallestAllowed = CoinMin(smallestAllowed,0.1*largestAllowed);
4692  randomNumberGenerator_.setSeed(1000000000);
4693#if ABC_NORMAL_DEBUG>0
4694  printf("perturbation %g constantPerturbation %g smallestNonZero %g maximumFraction %g smallestAllowed %g\n",
4695         perturbation,constantPerturbation,smallestNonZero,maximumFraction,
4696         smallestAllowed);
4697#endif
4698  for (int iSequence = maximumAbcNumberRows_; iSequence < numberTotal_; iSequence++) {
4699    if (lower[iSequence] < upper[iSequence] && getInternalStatus(iSequence) != basic) {
4700      double value = perturbation;
4701      double currentValue = cost[iSequence];
4702      double currentLargestAllowed=fabs(currentValue)*1.0e-3;
4703      if (!currentValue)
4704        currentLargestAllowed=1.0e-3;
4705      if (currentLargestAllowed<=10.0*smallestAllowed)
4706        continue; // don't perturb tiny value
4707      value = CoinMin(value, constantPerturbation + maximumFraction * (fabs(currentValue) + 1.0e-1 * perturbation + 1.0e-8));
4708      //value = CoinMin(value,constantPerturbation;+maximumFraction*fabs(currentValue));
4709      double value2 = constantPerturbation + 1.0e-1 * smallestNonZero;
4710      if (uniformChange) {
4711        value = maximumFraction;
4712        value2 = maximumFraction;
4713      }
4714      //double perturbationValue=2.0*(perturbationArray[iSequence]-0.5); // for now back to normal random
4715      double perturbationValue=randomNumberGenerator_.randomDouble();
4716      perturbationValue = 1.0-m1+m1*perturbationValue; // clean up
4717      if (lower[iSequence] > -largeValue_) {
4718        if (fabs(lower[iSequence]) <
4719            fabs(upper[iSequence])) {
4720          value *= perturbationValue;
4721          value2 *= perturbationValue;
4722        } else {
4723          //value *= -(1.0-m1+m1*randomNumberGenerator_.randomDouble());
4724          //value2 *= -(1.0-m1+m1*randomNumberGenerator_.randomDouble());
4725          value = 0.0;
4726        }
4727      } else if (upper[iSequence] < largeValue_) {
4728        value *= -perturbationValue;
4729        value2 *= -perturbationValue;
4730      } else {
4731        value = 0.0;
4732      }
4733      if (value) {
4734        int length = columnLength[iSequence];
4735        if (length > 3) {
4736          length = static_cast<int> (static_cast<double> (length) * factor);
4737          length = CoinMax(3, length);
4738        }
4739        double multiplier;
4740        if (length < 10)
4741          multiplier = weight[length];
4742        else
4743          multiplier = weight[10];
4744        value *= multiplier;
4745        value = CoinMin(value, value2);
4746        if (value) {
4747          // get in range
4748          if (fabs(value) <= smallestAllowed) {
4749            value *= 10.0;
4750            while (fabs(value) <= smallestAllowed)
4751              value *= 10.0;
4752          } else if (fabs(value) > currentLargestAllowed) {
4753            value *= 0.1;
4754            while (fabs(value) > currentLargestAllowed)
4755              value *= 0.1;
4756          }
4757        }
4758        if (currentValue) {
4759          largest = CoinMax(largest, fabs(value));
4760          if (fabs(value) > fabs(currentValue)*largestPerCent)
4761            largestPerCent = fabs(value / currentValue);
4762        } else {
4763          largestZero = CoinMax(largestZero, fabs(value));
4764        }
4765        // but negative if at ub
4766        if (getInternalStatus(iSequence) == atUpperBound)
4767          value = -value;
4768        if (printOut)
4769          printf("col %d cost %g change %g\n", iSequence, cost[iSequence], value);
4770        cost[iSequence] += value;
4771      }
4772    }
4773  }
4774#else
4775  // old way
4776     if (perturbation_ > 50) {
4777          // Experiment
4778          // maximumFraction could be 1.0e-10 to 1.0
4779          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};
4780          int whichOne = perturbation_ - 51;
4781          //if (inCbcOrOther&&whichOne>0)
4782          //whichOne--;
4783          maximumFraction = m[CoinMin(whichOne, 10)];
4784     } else if (inCbcOrOther) {
4785          //maximumFraction = 1.0e-6;
4786     }
4787     int iRow;
4788     double smallestNonZero = 1.0e100;
4789     numberNonZero = 0;
4790     if (perturbation_ >= 50) {
4791          perturbation = 1.0e-8;
4792          if (perturbation_ > 50 && perturbation_ < 60) 
4793            perturbation = CoinMax(1.0e-8,maximumFraction);
4794          bool allSame = true;
4795          double lastValue = 0.0;
4796          for (iRow = 0; iRow < numberRows_; iRow++) {
4797               double lo = abcLower_[iRow];
4798               double up = abcUpper_[iRow];
4799               if (lo && lo > -1.0e10) {
4800                    numberNonZero++;
4801                    lo = fabs(lo);
4802                    if (!lastValue)
4803                         lastValue = lo;
4804                    else if (fabs(lo - lastValue) > 1.0e-7)
4805                         allSame = false;
4806               }
4807               if (up && up < 1.0e10) {
4808                    numberNonZero++;
4809                    up = fabs(up);
4810                    if (!lastValue)
4811                         lastValue = up;
4812                    else if (fabs(up - lastValue) > 1.0e-7)
4813                         allSame = false;
4814               }
4815          }
4816          double lastValue2 = 0.0;
4817          for (int iColumn = maximumAbcNumberRows_; iColumn < numberTotal_; iColumn++) {
4818               double lo = abcLower_[iColumn];
4819               double up = abcUpper_[iColumn];
4820               if (lo < up) {
4821                    double value =
4822                         fabs(abcCost_[iColumn]);
4823                    perturbation = CoinMax(perturbation, value);
4824                    if (value) {
4825                         smallestNonZero = CoinMin(smallestNonZero, value);
4826                    }
4827               }
4828               if (lo && lo > -1.0e10) {
4829                    //numberNonZero++;
4830                    lo = fabs(lo);
4831                    if (!lastValue2)
4832                         lastValue2 = lo;
4833                    else if (fabs(lo - lastValue2) > 1.0e-7)
4834                         allSame = false;
4835               }
4836               if (up && up < 1.0e10) {
4837                    //numberNonZero++;
4838                    up = fabs(up);
4839                    if (!lastValue2)
4840                         lastValue2 = up;
4841                    else if (fabs(up - lastValue2) > 1.0e-7)
4842                         allSame = false;
4843               }
4844          }
4845          if (allSame) {
4846               // Check elements
4847               double smallestNegative;
4848               double largestNegative;
4849               double smallestPositive;
4850               double largestPositive;
4851               matrix_->rangeOfElements(smallestNegative, largestNegative,
4852                                        smallestPositive, largestPositive);
4853               if (smallestNegative == largestNegative &&
4854                         smallestPositive == largestPositive) {
4855                    // Really hit perturbation
4856                    double adjust = CoinMin(100.0 * maximumFraction, 1.0e-3 * CoinMax(lastValue, lastValue2));
4857                    maximumFraction = CoinMax(adjust, maximumFraction);
4858               }
4859          }
4860          perturbation = CoinMin(perturbation, smallestNonZero / maximumFraction);
4861     }
4862     double largestZero = 0.0;
4863     double largest = 0.0;
4864     double largestPerCent = 0.0;
4865     // modify costs
4866     bool printOut = (handler_->logLevel() == 63);
4867     printOut = false;
4868     // 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};
4869     // 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};
4870     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};
4871     //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};
4872     //double extraWeight = 10.0;
4873     // Scale back if wanted
4874     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};
4875     if (constantPerturbation < 99.0 * dualTolerance_) {
4876          perturbation *= 0.1;
4877          //extraWeight = 0.5;
4878          memcpy(weight, weight2, sizeof(weight2));
4879     }
4880     // adjust weights if all columns long
4881     double factor = 1.0;
4882     if (maxLength) {
4883          factor = 3.0 / static_cast<double> (minLength);
4884     }
4885     // Make variables with more elements more expensive
4886     const double m1 = 0.5;
4887     double smallestAllowed = CoinMin(1.0e-2 * dualTolerance_, maximumFraction);
4888     double largestAllowed = CoinMax(1.0e3 * dualTolerance_, maximumFraction * averageCost);
4889     // smaller if in BAB
4890     //if (inCbcOrOther)
4891     //largestAllowed=CoinMin(largestAllowed,1.0e-5);
4892     //smallestAllowed = CoinMin(smallestAllowed,0.1*largestAllowed);
4893     for (int iColumn = maximumAbcNumberRows_; iColumn < numberTotal_; iColumn++) {
4894          if (abcLower_[iColumn] < abcUpper_[iColumn] 
4895              && getInternalStatus(iColumn) != basic) {
4896               double value = perturbation;
4897               double currentValue = abcCost_[iColumn];
4898               value = CoinMin(value, constantPerturbation + maximumFraction * (fabs(currentValue) + 1.0e-1 * perturbation + 1.0e-8));
4899               //value = CoinMin(value,constantPerturbation;+maximumFraction*fabs(currentValue));
4900               double value2 = constantPerturbation + 1.0e-1 * smallestNonZero;
4901               if (uniformChange) {
4902                    value = maximumFraction;
4903                    value2 = maximumFraction;
4904               }
4905               if (abcLower_[iColumn] > -largeValue_) {
4906                    if (fabs(abcLower_[iColumn]) <
4907                              fabs(abcUpper_[iColumn])) {
4908                         value *= (1.0 - m1 + m1 * randomNumberGenerator_.randomDouble());
4909                         value2 *= (1.0 - m1 + m1 * randomNumberGenerator_.randomDouble());
4910                    } else {
4911                         //value *= -(1.0-m1+m1*randomNumberGenerator_.randomDouble());
4912                         //value2 *= -(1.0-m1+m1*randomNumberGenerator_.randomDouble());
4913                         value = 0.0;
4914                    }
4915               } else if (abcUpper_[iColumn] < largeValue_) {
4916                    value *= -(1.0 - m1 + m1 * randomNumberGenerator_.randomDouble());
4917                    value2 *= -(1.0 - m1 + m1 * randomNumberGenerator_.randomDouble());
4918               } else {
4919                    value = 0.0;
4920               }
4921               if (value) {
4922                    int length = columnLength[iColumn];
4923                    if (length > 3) {
4924                         length = static_cast<int> (static_cast<double> (length) * factor);
4925                         length = CoinMax(3, length);
4926                    }
4927                    double multiplier;
4928#if 1
4929                    if (length < 10)
4930                         multiplier = weight[length];
4931                    else
4932                         multiplier = weight[10];
4933#else
4934                    if (length < 10)
4935                         multiplier = weight[length];
4936                    else
4937                         multiplier = weight[10] + extraWeight * (length - 10);
4938                    multiplier *= 0.5;
4939#endif
4940                    value *= multiplier;
4941                    value = CoinMin(value, value2);
4942                    if (value) {
4943                         // get in range
4944                         if (fabs(value) <= smallestAllowed) {
4945                              value *= 10.0;
4946                              while (fabs(value) <= smallestAllowed)
4947                                   value *= 10.0;
4948                         } else if (fabs(value) > largestAllowed) {
4949                              value *= 0.1;
4950                              while (fabs(value) > largestAllowed)
4951                                   value *= 0.1;
4952                         }
4953                    }
4954                    if (currentValue) {
4955                         largest = CoinMax(largest, fabs(value));
4956                         if (fabs(value) > fabs(currentValue)*largestPerCent)
4957                              largestPerCent = fabs(value / currentValue);
4958                    } else {
4959                         largestZero = CoinMax(largestZero, fabs(value));
4960                    }
4961                    // but negative if at ub
4962                    if (getInternalStatus(iColumn) == atUpperBound)
4963                         value = -value;
4964                    if (printOut)
4965                         printf("col %d cost %g change %g\n", iColumn, objectiveWork_[iColumn], value);
4966                    abcCost_[iColumn] += value;
4967               }
4968          }
4969     }
4970#endif
4971#ifdef PERT_STATISTICS
4972  {
4973    double averageCost = 0.0;
4974    int numberNonZero = 0;
4975    double * COIN_RESTRICT sort = new double[numberColumns_];
4976    for (int i = 0; i < numberColumns_; i++) {
4977      double value = fabs(abcCost_[i+maximumAbcNumberRows_]);
4978      sort[i] = value;
4979      averageCost += value;
4980      if (value)
4981        numberNonZero++;
4982    }
4983    if (numberNonZero)
4984      averageCost /= static_cast<double> (numberNonZero);
4985    else
4986      averageCost = 1.0;
4987    std::sort(sort, sort + numberColumns_);
4988    int number = 1;
4989    double last = sort[0];
4990    for (int i = 1; i < numberColumns_; i++) {
4991      if (last != sort[i])
4992        number++;
4993    last = sort[i];
4994    }
4995    delete [] sort;
4996#if ABC_NORMAL_DEBUG>0
4997    printf("nnz %d percent %d", number, (number * 100) / numberColumns_);
4998    breakdown("Objective after perturbation", numberColumns_, abcCost_+numberRows_);
4999#endif
5000  }
5001#endif
5002  perturbation_=100;
5003  // copy cost to  perturbation
5004#ifdef CLEANUP_DJS
5005  CoinAbcMemcpy(abcPerturbation_,abcCost_,numberTotal_);
5006#endif
5007  handler_->message(CLP_SIMPLEX_PERTURB, messages_)
5008    << 100.0 * maximumFraction << perturbation << largest << 100.0 * largestPerCent << largestZero
5009    << CoinMessageEol;
5010}
5011/* Does something about fake tolerances
5012   -1 initializes
5013   1 might be optimal (back to original costs)
5014   2 after change currentDualBound
5015   3 might be optimal BUT change currentDualBound
5016*/
5017int 
5018AbcSimplexDual::bounceTolerances(int type)
5019{
5020  // for now turn off
5021  if (type==-1) {
5022    if (!perturbationFactor_)
5023      perturbationFactor_=0.001;
5024    CoinAbcMemcpy(abcCost_,costSaved_,numberTotal_);
5025    if (!perturbationFactor_) {
5026      //currentDualBound_=CoinMin(dualBound_,1.0e11);
5027      //CoinAbcMultiplyAdd(perturbationSaved_,numberTotal_,currentDualTolerance_,
5028      //                 abcPerturbation_,0.0);
5029    } else {
5030      //currentDualBound_=CoinMin(dualBound_,1.0e11);
5031      double factor=currentDualTolerance_*perturbationFactor_;
5032      perturbB(factor,0);
5033    }
5034#if ABC_NORMAL_DEBUG>0
5035    //printf("new dual bound %g (bounceTolerances)\n",currentDualBound_);
5036#endif
5037  } else {
5038    if (stateDualColumn_==-2) {
5039      double newTolerance=CoinMin(1.0e-1,dualTolerance_*100.0);
5040#if ABC_NORMAL_DEBUG>0
5041      printf("After %d iterations reducing dual tolerance from %g to %g\n",
5042             numberIterations_,currentDualTolerance_,newTolerance);
5043#endif
5044      currentDualTolerance_=newTolerance;
5045      stateDualColumn_=-1;
5046    } else if (stateDualColumn_==-1) {
5047#if ABC_NORMAL_DEBUG>0
5048      printf("changing dual tolerance from %g to %g (vaguely optimal)\n",
5049             currentDualTolerance_,dualTolerance_);
5050#endif
5051      currentDualTolerance_=dualTolerance_;
5052      stateDualColumn_=0;
5053    } else {
5054      if (stateDualColumn_>1) {
5055        normalDualColumnIteration_ = numberIterations_+CoinMax(10,numberRows_/4);
5056#if ABC_NORMAL_DEBUG>0
5057        printf("no cost modification until iteration %d\n",normalDualColumnIteration_);
5058#endif
5059      }
5060    }
5061    resetFakeBounds(0);
5062    CoinAbcMemcpy(abcCost_,abcPerturbation_,numberTotal_);
5063    numberChanged_ = 0; // Number of variables with changed costs
5064    moveToBasic();
5065    if (type!=1)
5066      return 0;
5067    // save status and solution
5068    CoinAbcMemcpy(status_,internalStatus_,numberTotal_);
5069    CoinAbcMemcpy(solutionSaved_,abcSolution_,numberTotal_);
5070#if ABC_NORMAL_DEBUG>0
5071    printf("A sumDual %g sumPrimal %g\n",sumDualInfeasibilities_,sumPrimalInfeasibilities_);
5072#endif
5073    gutsOfSolution(3);
5074#if ABC_NORMAL_DEBUG>0
5075    printf("B sumDual %g sumPrimal %g\n",sumDualInfeasibilities_,sumPrimalInfeasibilities_);
5076#endif
5077    double saveSum1=sumPrimalInfeasibilities_;
5078    double saveSumD1=sumDualInfeasibilities_;
5079    bool goToPrimal=false;
5080    if (!sumOfRelaxedDualInfeasibilities_) {
5081      if (perturbation_<101) {
5082        // try going back to original costs
5083        CoinAbcMemcpy(abcCost_,costSaved_,numberTotal_);
5084        perturbation_=101;
5085        numberChanged_ = 0; // Number of variables with changed costs
5086        moveToBasic();
5087        gutsOfSolution(3);
5088#if ABC_NORMAL_DEBUG>0
5089        printf("B2 sumDual %g sumPrimal %g\n",sumDualInfeasibilities_,sumPrimalInfeasibilities_);
5090#endif
5091        saveSumD1=sumDualInfeasibilities_;
5092      }
5093      if (sumOfRelaxedDualInfeasibilities_) {
5094        makeNonFreeVariablesDualFeasible();
5095        moveToBasic();
5096        abcProgress_.reset();
5097        gutsOfSolution(3);
5098#if ABC_NORMAL_DEBUG>0
5099        printf("C sumDual %g sumPrimal %g\n",sumDualInfeasibilities_,sumPrimalInfeasibilities_);
5100#endif
5101        if (sumPrimalInfeasibilities_>1.0e2*saveSum1+currentDualBound_*0.2
5102          /*|| sumOfRelaxedDualInfeasibilities_<1.0e-4*sumPrimalInfeasibilities_*/) {
5103          numberDisasters_++;
5104          if (numberDisasters_>2)
5105            goToPrimal=true;
5106        }
5107      } else {
5108        numberDualInfeasibilities_=0; 
5109        sumDualInfeasibilities_=0.0; 
5110      }
5111    } else {
5112      // carry on
5113      makeNonFreeVariablesDualFeasible();
5114      moveToBasic();
5115      abcProgress_.reset();
5116      gutsOfSolution(3);
5117#if ABC_NORMAL_DEBUG>0
5118      printf("C2 sumDual %g sumPrimal %g\n",sumDualInfeasibilities_,sumPrimalInfeasibilities_);
5119#endif
5120      numberDisasters_++;
5121      if (sumPrimalInfeasibilities_>1.0e2*saveSum1+currentDualBound_*0.2
5122          ||(saveSumD1<1.0e-4*sumPrimalInfeasibilities_&&saveSumD1<1.0)) {
5123        numberDisasters_++;
5124        if (numberDisasters_>2||saveSumD1<1.0e-10*sumPrimalInfeasibilities_)
5125          goToPrimal=true;
5126      }
5127    }
5128    if (goToPrimal) {
5129      // go to primal
5130      CoinAbcMemcpy(internalStatus_,status_,numberTotal_);
5131      CoinAbcMemcpy(abcSolution_,solutionSaved_,numberTotal_);
5132      int *  COIN_RESTRICT abcPivotVariable = abcPivotVariable_;
5133      // redo pivotVariable
5134      int numberBasic=0;
5135      for (int i=0;i<numberTotal_;i++) {
5136        if (getInternalStatus(i)==basic)
5137          abcPivotVariable[numberBasic++]=i;
5138      }
5139      assert (numberBasic==numberRows_);
5140      moveToBasic();
5141      checkPrimalSolution(false);
5142#if ABC_NORMAL_DEBUG>0
5143      printf("D sumPrimal %g\n",sumPrimalInfeasibilities_);
5144#endif
5145      return 1;
5146    }
5147  }
5148
5149  CoinAbcGatherFrom(abcCost_,costBasic_,abcPivotVariable_,numberRows_);
5150  return 0;
5151}
5152#if ABC_NORMAL_DEBUG>0
5153extern int cilk_conflict;
5154#endif
5155int 
5156AbcSimplexDual::noPivotRow()
5157{
5158#if ABC_NORMAL_DEBUG>3
5159  if (handler_->logLevel() & 32)
5160    printf("** no row pivot\n");
5161#endif
5162  int numberPivots = abcFactorization_->pivots();
5163  bool specialCase;
5164  int useNumberFake;
5165  int returnCode = 0;
5166  if (numberPivots <= CoinMax(dontFactorizePivots_, 20) &&
5167      (specialOptions_ & 2048) != 0 && currentDualBound_ >= 1.0e8) {
5168    specialCase = true;
5169    // as dual bound high - should be okay
5170    useNumberFake = 0;
5171  } else {
5172    specialCase = false;
5173    useNumberFake = numberAtFakeBound();
5174  }
5175  if (!numberPivots || specialCase) {
5176    // may have crept through - so may be optimal
5177    // check any flagged variables
5178    if (numberFlagged_ && numberPivots) {
5179      // try factorization
5180      returnCode = -2;
5181    }
5182   
5183    if (useNumberFake || numberDualInfeasibilities_) {
5184      // may be dual infeasible
5185      if ((specialOptions_ & 1024) == 0)
5186        problemStatus_ = -5;
5187      else if (!useNumberFake && numberPrimalInfeasibilities_
5188               && !numberPivots)
5189        problemStatus_ = 1;
5190    } else {
5191      if (numberFlagged_) {
5192#if ABC_NORMAL_DEBUG>3
5193        std::cout << "Flagged variables at end - infeasible?" << std::endl;
5194        printf("Probably infeasible - pivot was %g\n", alpha_);
5195#endif
5196#if ABC_NORMAL_DEBUG>3
5197        abort();
5198#endif
5199        problemStatus_ = -5;
5200      } else {
5201        problemStatus_ = 0;
5202#ifndef CLP_CHECK_NUMBER_PIVOTS
5203#define CLP_CHECK_NUMBER_PIVOTS 10
5204#endif
5205#if CLP_CHECK_NUMBER_PIVOTS < 20
5206        if (numberPivots > CLP_CHECK_NUMBER_PIVOTS) {
5207          gutsOfSolution(2);
5208          rawObjectiveValue_ +=sumNonBasicCosts_;
5209          setClpSimplexObjectiveValue();
5210          if (numberPrimalInfeasibilities_) {
5211            problemStatus_ = -1;
5212            returnCode = -2;
5213          }
5214        }
5215#endif
5216        if (!problemStatus_) {
5217          // make it look OK
5218          numberPrimalInfeasibilities_ = 0;
5219          sumPrimalInfeasibilities_ = 0.0;
5220          numberDualInfeasibilities_ = 0;
5221          sumDualInfeasibilities_ = 0.0;
5222          if (numberChanged_) {
5223            numberChanged_ = 0; // Number of variables with changed costs
5224            CoinAbcMemcpy(abcCost_,abcCost_+maximumNumberTotal_,numberTotal_);
5225            CoinAbcGatherFrom(abcCost_,costBasic_,abcPivotVariable_,numberRows_);
5226                  // make sure duals are current
5227            gutsOfSolution(1);
5228            setClpSimplexObjectiveValue();
5229            abcProgress_.modifyObjective(-COIN_DBL_MAX);
5230            if (numberDualInfeasibilities_) {
5231              problemStatus_ = 10; // was -3;
5232              numberTimesOptimal_++;
5233            } else {
5234              computeObjectiveValue(true);
5235            }
5236          } else if (numberPivots) {
5237            computeObjectiveValue(true);
5238          }
5239          sumPrimalInfeasibilities_ = 0.0;
5240        }
5241        if ((specialOptions_&(1024 + 16384)) != 0 && !problemStatus_&&abcFactorization_->pivots()) {
5242          CoinIndexedVector * arrayVector = &usefulArray_[arrayForFtran_];
5243          double * rhs = arrayVector->denseVector();
5244          arrayVector->clear();
5245          CoinAbcScatterTo(solutionBasic_,abcSolution_,abcPivotVariable_,numberRows_);
5246          abcMatrix_->timesModifyExcludingSlacks(-1.0, abcSolution_,rhs);
5247          //#define CHECK_ACCURACY
5248#ifdef CHECK_ACCURACY
5249          bool bad = false;
5250#endif
5251          bool bad2 = false;
5252          int i;
5253          for ( i = 0; i < numberRows_; i++) {
5254            if (rhs[i] < abcLower_[i] - primalTolerance_ ||
5255                rhs[i] > abcUpper_[i] + primalTolerance_) {
5256              bad2 = true;
5257#ifdef CHECK_ACCURACY
5258              printf("row %d out of bounds %g, %g correct %g bad %g\n", i,
5259                     abcLower_[i], abcUpper_[i],
5260                     rhs[i], abcSolution_[i]);
5261#endif
5262            } else if (fabs(rhs[i] - abcSolution_[i]) > 1.0e-3) {
5263#ifdef CHECK_ACCURACY
5264              bad = true;
5265              printf("row %d correct %g bad %g\n", i, rhs[i], abcSolution_[i]);
5266#endif
5267            }
5268            rhs[i] = 0.0;
5269          }
5270          for (int i = 0; i < numberTotal_; i++) {
5271            if (abcSolution_[i] < abcLower_[i] - primalTolerance_ ||
5272                abcSolution_[i] > abcUpper_[i] + primalTolerance_) {
5273              bad2 = true;
5274#ifdef CHECK_ACCURACY
5275              printf("column %d out of bounds %g, %g value %g\n", i,
5276                     abcLower_[i], abcUpper_[i],
5277                     abcSolution_[i]);
5278#endif
5279            }
5280          }
5281          if (bad2) {
5282            problemStatus_ = -3;
5283            returnCode = -2;
5284            // Force to re-factorize early next time
5285            int numberPivots = abcFactorization_->pivots();
5286            if (forceFactorization_<0)
5287              forceFactorization_=100000;
5288            forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
5289          }
5290        }
5291      }
5292    }
5293  } else {
5294    problemStatus_ = -3;
5295    returnCode = -2;
5296    // Force to re-factorize early next time
5297    int numberPivots = abcFactorization_->pivots();
5298    if (forceFactorization_<0)
5299      forceFactorization_=100000;
5300    forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
5301  }
5302  return returnCode;
5303}
5304void 
5305AbcSimplexDual::dualPivotColumn()
5306{
5307  double saveAcceptable=acceptablePivot_;
5308  if (largestPrimalError_>1.0e-5||largestDualError_>1.0e-5) {
5309    //if ((largestPrimalError_>1.0e-5||largestDualError_>1.0e-5)&&false) {
5310    if (!abcFactorization_->pivots())
5311      acceptablePivot_*=1.0e2;
5312    else if (abcFactorization_->pivots()<5)
5313      acceptablePivot_*=1.0e1;
5314  }
5315  dualColumn2();
5316  if (sequenceIn_<0&&acceptablePivot_>saveAcceptable) {
5317    acceptablePivot_=saveAcceptable;
5318#if ABC_NORMAL_DEBUG>0
5319    printf("No pivot!!!\n");
5320#endif
5321    // try again
5322    dualColumn1(true);
5323    dualColumn2();
5324  } else {
5325    acceptablePivot_=saveAcceptable;
5326  }
5327  if ((stateOfProblem_&VALUES_PASS)!=0) {
5328    // check no flips
5329    assert (!usefulArray_[arrayForFlipBounds_].getNumElements());
5330    // If no pivot - then try other way if plausible
5331    if (sequenceIn_<0) {
5332      // could use fake primal basic values
5333      if ((directionOut_<0&&lowerOut_>-1.0e30)||
5334          (directionOut_>0&&upperOut_<1.30)) {
5335        directionOut_=-directionOut_;
5336        if (directionOut_<0&&abcDj_[sequenceOut_]<0.0)
5337          upperTheta_=-abcDj_[sequenceOut_];
5338        else if (directionOut_>0&&abcDj_[sequenceOut_]>0.0)
5339          upperTheta_=abcDj_[sequenceOut_];
5340        else
5341          upperTheta_=1.0e30;
5342        CoinPartitionedVector & tableauRow = usefulArray_[arrayForTableauRow_];
5343        CoinPartitionedVector & candidateList = usefulArray_[arrayForDualColumn_];
5344        assert (!candidateList.getNumElements());
5345        tableauRow.compact();
5346        candidateList.compact();
5347        // redo candidate list
5348        int *  COIN_RESTRICT index = tableauRow.getIndices();
5349        double *  COIN_RESTRICT array = tableauRow.denseVector();
5350        double *  COIN_RESTRICT arrayCandidate=candidateList.denseVector();
5351        int *  COIN_RESTRICT indexCandidate = candidateList.getIndices();
5352        const double *  COIN_RESTRICT abcDj = abcDj_;
5353        const unsigned char *  COIN_RESTRICT internalStatus = internalStatus_;
5354        // do first pass to get possibles
5355        double bestPossible = 0.0;
5356        // We can also see if infeasible or pivoting on free
5357        double tentativeTheta = 1.0e25; // try with this much smaller as guess
5358        double acceptablePivot = currentAcceptablePivot_;
5359        double dualT=-currentDualTolerance_;
5360        // fixed will have been taken out by now
5361        const double multiplier[] = { 1.0, -1.0};
5362        int freeSequence=-1;
5363        double freeAlpha=0.0;
5364        int numberNonZero=tableauRow.getNumElements();
5365        int numberRemaining=0;
5366        if (ordinaryVariables()) {
5367          for (int i = 0; i < numberNonZero; i++) {
5368            int iSequence = index[i];
5369            double tableauValue=array[i];
5370            unsigned char iStatus=internalStatus[iSequence]&7;
5371            double mult = multiplier[iStatus];
5372            double alpha = tableauValue * mult;
5373            double oldValue = abcDj[iSequence] * mult;
5374            double value = oldValue - tentativeTheta * alpha;
5375            if (value < dualT) {
5376              bestPossible = CoinMax(bestPossible, alpha);
5377              value = oldValue - upperTheta_ * alpha;
5378              if (value < dualT && alpha >= acceptablePivot) {
5379                upperTheta_ = (oldValue - dualT) / alpha;
5380              }
5381              // add to list
5382              arrayCandidate[numberRemaining] = alpha;
5383              indexCandidate[numberRemaining++] = iSequence;
5384            }
5385          }
5386        } else {
5387          double badFree = 0.0;
5388          freeAlpha = currentAcceptablePivot_;
5389          double currentDualTolerance = currentDualTolerance_;
5390          for (int i = 0; i < numberNonZero; i++) {
5391            int iSequence = index[i];
5392            double tableauValue=array[i];
5393            unsigned char iStatus=internalStatus[iSequence]&7;
5394            if ((iStatus&4)==0) {
5395              double mult = multiplier[iStatus];
5396              double alpha = tableauValue * mult;
5397              double oldValue = abcDj[iSequence] * mult;
5398              double value = oldValue - tentativeTheta * alpha;
5399              if (value < dualT) {
5400                bestPossible = CoinMax(bestPossible, alpha);
5401                value = oldValue - upperTheta_ * alpha;
5402                if (value < dualT && alpha >= acceptablePivot) {
5403                  upperTheta_ = (oldValue - dualT) / alpha;
5404                }
5405                // add to list
5406                arrayCandidate[numberRemaining] = alpha;
5407                indexCandidate[numberRemaining++] = iSequence;
5408              }
5409            } else {
5410              bool keep;
5411              bestPossible = CoinMax(bestPossible, fabs(tableauValue));
5412              double oldValue = abcDj[iSequence];
5413              // If free has to be very large - should come in via dualRow
5414              //if (getInternalStatus(iSequence+addSequence)==isFree&&fabs(tableauValue)<1.0e-3)
5415              //break;
5416              if (oldValue > currentDualTolerance) {
5417                keep = true;
5418              } else if (oldValue < -currentDualTolerance) {
5419                keep = true;
5420              } else {
5421                if (fabs(tableauValue) > CoinMax(10.0 * acceptablePivot, 1.0e-5)) {
5422                  keep = true;
5423                } else {
5424                  keep = false;
5425                  badFree = CoinMax(badFree, fabs(tableauValue));
5426                }
5427              }
5428              if (keep) {
5429                // free - choose largest
5430                if (fabs(tableauValue) > fabs(freeAlpha)) {
5431                  freeAlpha = tableauValue;
5432                  freeSequence = iSequence;
5433                }
5434              }
5435            }
5436          }
5437        }
5438        candidateList.setNumElements(numberRemaining);
5439        if (freeSequence>=0) {
5440          sequenceIn_=freeSequence;
5441        } else if (fabs(upperTheta_-fabs(abcDj_[sequenceOut_]))<dualTolerance_) {
5442          // can just move
5443          stateOfIteration_=-1;
5444          return;
5445        }
5446        dualColumn2();
5447      }
5448    }
5449    // redo dualOut_
5450    if (directionOut_<0) {
5451      dualOut_ = valueOut_ - upperOut_;
5452    } else {
5453      dualOut_ = lowerOut_ - valueOut_;
5454    }
5455  }
5456  usefulArray_[arrayForDualColumn_].clear();
5457  //clearArrays(arrayForDualColumn_);
5458  if (sequenceIn_ >= 0&&upperTheta_==COIN_DBL_MAX) {
5459    lowerIn_ = abcLower_[sequenceIn_];
5460    upperIn_ = abcUpper_[sequenceIn_];
5461    valueIn_ = abcSolution_[sequenceIn_];
5462    dualIn_ = abcDj_[sequenceIn_];
5463    if (valueIn_<lowerIn_+primalTolerance_) {
5464      directionIn_ = 1;
5465    } else if (valueIn_>upperIn_-primalTolerance_) {
5466      directionIn_ = -1;
5467    } else {
5468      if (alpha_ < 0.0) {
5469        // as if from upper bound
5470        directionIn_ = -1;
5471        //assert(upperIn_ == valueIn_);
5472      } else {
5473        // as if from lower bound
5474        directionIn_ = 1;
5475        //assert(lowerIn_ == valueIn_);
5476      }
5477    }
5478  }
5479#if ABC_DEBUG
5480  checkArrays(4);
5481#endif
5482#if CAN_HAVE_ZERO_OBJ>1
5483  if ((specialOptions_&2097152)!=0) 
5484    theta_=0.0;
5485#endif
5486  movement_=0.0;
5487  objectiveChange_=0.0;
5488  /*
5489    0 - take iteration
5490    1 - don't take but continue
5491    2 - don't take and break
5492  */
5493  btranAlpha_ = -alpha_ * directionOut_; // for check
5494}
5495void
5496AbcSimplexDual::getTableauColumnPart2()
5497{
5498#if ABC_PARALLEL
5499  if (parallelMode_!=0) {
5500    abcFactorization_->updateColumnFTPart2(usefulArray_[arrayForFtran_]);
5501    // pivot element
5502    alpha_ = usefulArray_[arrayForFtran_].denseVector()[pivotRow_];
5503  }
5504#endif
5505}
5506void 
5507AbcSimplexDual::replaceColumnPart3()
5508{
5509  abcFactorization_->replaceColumnPart3(this,
5510                                        &usefulArray_[arrayForReplaceColumn_],
5511                                        &usefulArray_[arrayForFtran_],
5512                                        lastPivotRow_,
5513                                        ftAlpha_?ftAlpha_:alpha_);
5514}
5515void
5516AbcSimplexDual::checkReplacePart1()
5517{
5518  //abcFactorization_->checkReplacePart1a(&usefulArray_[arrayForReplaceColumn_],pivotRow_);
5519  //usefulArray_[arrayForReplaceColumn_].print();
5520  ftAlpha_=abcFactorization_->checkReplacePart1(&usefulArray_[arrayForReplaceColumn_],pivotRow_);
5521}
5522#if 0
5523void
5524AbcSimplexDual::checkReplacePart1a()
5525{
5526  abcFactorization_->checkReplacePart1a(&usefulArray_[arrayForReplaceColumn_],pivotRow_);
5527}
5528void
5529AbcSimplexDual::checkReplacePart1b()
5530{
5531  //usefulArray_[arrayForReplaceColumn_].print();
5532  ftAlpha_=abcFactorization_->checkReplacePart1b(&usefulArray_[arrayForReplaceColumn_],pivotRow_);
5533}
5534#endif
5535/*
5536  returns
5537  0 - OK
5538  1 - take iteration then re-factorize
5539  2 - flag something and skip
5540  3 - break and re-factorize
5541  5 - take iteration then re-factorize because of memory
5542 */
5543int 
5544AbcSimplexDual::checkReplace()
5545{
5546  int returnCode=0;
5547  if (!stateOfIteration_) {
5548    int saveStatus=problemStatus_;
5549    // check update
5550#if MOVE_REPLACE_PART1A < 0
5551    checkReplacePart1();
5552#endif
5553    int updateStatus=abcFactorization_->checkReplacePart2(pivotRow_,btranAlpha_,alpha_,ftAlpha_);
5554#if ABC_DEBUG
5555    checkArrays();
5556#endif
5557    // If looks like bad pivot - refactorize
5558    if (fabs(dualOut_) > 1.0e50)
5559      updateStatus = 2;
5560    // if no pivots, bad update but reasonable alpha - take and invert
5561    if (updateStatus == 2 &&
5562        !abcFactorization_->pivots() && fabs(alpha_) > 1.0e-5)
5563      updateStatus = 4;
5564    if (updateStatus == 1 || updateStatus == 4) {
5565      // slight error
5566      if (abcFactorization_->pivots() > 5 || updateStatus == 4) {
5567        problemStatus_ = -2; // factorize now
5568        returnCode = -3;
5569      }
5570    } else if (updateStatus == 2) {
5571      // major error
5572      // later we may need to unwind more e.g. fake bounds
5573      if (abcFactorization_->pivots() &&
5574          ((moreSpecialOptions_ & 16) == 0 || abcFactorization_->pivots() > 4)) {
5575        problemStatus_ = -2; // factorize now
5576        returnCode = -2;
5577        moreSpecialOptions_ |= 16;
5578        stateOfIteration_=2;
5579      } else {
5580        assert (sequenceOut_>=0);
5581        // need to reject something
5582        char x = isColumn(sequenceOut_) ? 'C' : 'R';
5583        handler_->message(CLP_SIMPLEX_FLAG, messages_)
5584          << x << sequenceWithin(sequenceOut_)
5585          << CoinMessageEol;
5586#if ABC_NORMAL_DEBUG>0
5587        printf("BAD flag b %g pivotRow %d (seqout %d) sequenceIn %d\n", alpha_,pivotRow_,
5588               sequenceOut_,sequenceIn_);
5589#endif
5590        setFlagged(sequenceOut_);
5591        // so can't happen immediately again
5592        sequenceOut_=-1;
5593        //abcProgress_.clearBadTimes();
5594        lastBadIteration_ = numberIterations_; // say be more cautious
5595        //clearArrays(-1);
5596        //if(primalSolutionUpdated)
5597        //gutsOfSolution(2);
5598        stateOfIteration_=1;
5599      }
5600    } else if (updateStatus == 3) {
5601      // out of memory
5602      // increase space if not many iterations
5603