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

Last change on this file since 1910 was 1910, checked in by stefan, 7 years ago
  • add configure option --enable-aboca={1,2,3,4,yes,no}
  • compile Aboca source only if --enable-aboca set (instead of compiling empty source files)
  • fix svn properties
  • Property svn:keywords set to Id
File size: 201.6 KB
Line 
1/* $Id: AbcSimplexDual.cpp 1910 2013-01-27 02:00:13Z stefan $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others, Copyright (C) 2012, FasterCoin.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6/* Notes on implementation of dual simplex algorithm.
7   
8   When dual feasible:
9   
10   If primal feasible, we are optimal.  Otherwise choose an infeasible
11   basic variable to leave basis (normally going to nearest bound) (B).  We
12   now need to find an incoming variable which will leave problem
13   dual feasible so we get the row of the tableau corresponding to
14   the basic variable (with the correct sign depending if basic variable
15   above or below feasibility region - as that affects whether reduced
16   cost on outgoing variable has to be positive or negative).
17   
18   We now perform a ratio test to determine which incoming variable will
19   preserve dual feasibility (C).  If no variable found then problem
20   is infeasible (in primal sense).  If there is a variable, we then
21   perform pivot and repeat.  Trivial?
22   
23   -------------------------------------------
24   
25   A) How do we get dual feasible?  If all variables have bounds then
26   it is trivial to get feasible by putting non-basic variables to
27   correct bounds.  OSL did not have a phase 1/phase 2 approach but
28   instead effectively put fake bounds on variables and this is the
29   approach here, although I had hoped to make it cleaner.
30   
31   If there is a weight of X on getting dual feasible:
32   Non-basic variables with negative reduced costs are put to
33   lesser of their upper bound and their lower bound + X.
34   Similarly, mutatis mutandis, for positive reduced costs.
35   
36   Free variables should normally be in basis, otherwise I have
37   coding which may be able to come out (and may not be correct).
38   
39   In OSL, this weight was changed heuristically, here at present
40   it is only increased if problem looks finished.  If problem is
41   feasible I check for unboundedness.  If not unbounded we
42   could play with going into primal.  As long as weights increase
43   any algorithm would be finite.
44   
45   B) Which outgoing variable to choose is a virtual base class.
46   For difficult problems steepest edge is preferred while for
47   very easy (large) problems we will need partial scan.
48   
49   C) Sounds easy, but this is hardest part of algorithm.
50   1) Instead of stopping at first choice, we may be able
51   to flip that variable to other bound and if objective
52   still improving choose again.  These mini iterations can
53   increase speed by orders of magnitude but we may need to
54   go to more of a bucket choice of variable rather than looking
55   at them one by one (for speed).
56   2) Accuracy.  Reduced costs may be of wrong sign but less than
57   tolerance.  Pivoting on these makes objective go backwards.
58   OSL modified cost so a zero move was made, Gill et al
59   (in primal analogue) modified so a strictly positive move was
60   made.  It is not quite as neat in dual but that is what we
61   try and do.  The two problems are that re-factorizations can
62   change reduced costs above and below tolerances and that when
63   finished we need to reset costs and try again.
64   3) Degeneracy.  Gill et al helps but may not be enough.  We
65   may need more.  Also it can improve speed a lot if we perturb
66   the costs significantly.
67   
68   References:
69   Forrest and Goldfarb, Steepest-edge simplex algorithms for
70   linear programming - Mathematical Programming 1992
71   Forrest and Tomlin, Implementing the simplex method for
72   the Optimization Subroutine Library - IBM Systems Journal 1992
73   Gill, Murray, Saunders, Wright A Practical Anti-Cycling
74   Procedure for Linear and Nonlinear Programming SOL report 1988
75   
76   
77   TODO:
78   
79   a) Better recovery procedures.  At present I never check on forward
80   progress.  There is checkpoint/restart with reducing
81   re-factorization frequency, but this is only on singular
82   factorizations.
83   b) Fast methods for large easy problems (and also the option for
84   the code to automatically choose which method).
85   c) We need to be able to stop in various ways for OSI - this
86   is fairly easy.
87   
88*/
89
90#include "CoinPragma.hpp"
91
92#include <math.h>
93
94#include "CoinHelperFunctions.hpp"
95#include "CoinAbcHelperFunctions.hpp"
96#include "AbcSimplexDual.hpp"
97#include "ClpEventHandler.hpp"
98#include "AbcSimplexFactorization.hpp"
99#include "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 lastPivotValue = 0.0;
1598  double thisPivotValue=0.0;
1599  // if free then always choose , otherwise ...
1600  double theta = 1.0e50;
1601  // now flip flop between spare arrays until reasonable theta
1602  tentativeTheta = CoinMax(10.0 * upperTheta_, 1.0e-7);
1603 
1604  // loops increasing tentative theta until can't go through
1605#ifdef PRINT_RATIO_PROGRESS
1606  int iPass=0;
1607#endif
1608  dualColumnResult result2[8];
1609  for (int iBlock=0;iBlock<numberBlocks;iBlock++) {
1610    result2[iBlock].block=iBlock;
1611    if (numberBlocksReal)
1612      result2[iBlock].numberRemaining=candidateList.getNumElements(iBlock);
1613    else
1614      result2[iBlock].numberRemaining=candidateList.getNumElements();
1615    result2[iBlock].numberSwapped=0;
1616    result2[iBlock].numberLastSwapped=0;
1617  }
1618  double useThru=0.0;
1619  while (tentativeTheta < 1.0e22) {
1620    double thruThis = 0.0;
1621   
1622    double bestPivot = currentAcceptablePivot_;
1623    int bestSequence = -1;
1624    sequenceIn=-1;
1625    thisPivotValue=0.0;
1626    int numberInList=numberRemaining;
1627    double upperTheta = 1.0e50;
1628    double increaseInThis = 0.0; //objective increase in this loop
1629    for (int iBlock=0;iBlock<numberBlocks;iBlock++) {
1630      result2[iBlock].tentativeTheta=tentativeTheta;
1631      //result2[iBlock].numberRemaining=numberRemaining;
1632      //result2[iBlock].numberSwapped=numberSwapped;
1633      //result2[iBlock].numberLastSwapped=numberLastSwapped;
1634    }
1635    for (int iBlock=1;iBlock<numberBlocks;iBlock++) {
1636      cilk_spawn dualColumn2First(result2[iBlock]);
1637    }
1638    dualColumn2First(result2[0]);
1639    cilk_sync;
1640    //dualColumn2Bit(this,result2,numberBlocks);
1641    numberSwapped=0;
1642    numberRemaining=0;
1643    for (int iBlock=0;iBlock<numberBlocks;iBlock++) {
1644      if (result2[iBlock].bestEverPivot>bestPivot) {
1645        bestPivot=result2[iBlock].bestEverPivot;
1646        bestSequence=result2[iBlock].sequence;
1647      }
1648      numberSwapped+=result2[iBlock].numberSwapped;
1649      numberRemaining+=result2[iBlock].numberRemaining;
1650      thruThis+=result2[iBlock].thruThis;
1651      increaseInThis+=result2[iBlock].increaseInThis;
1652      if (upperTheta>result2[iBlock].theta)
1653        upperTheta=result2[iBlock].theta;
1654    }
1655    // end of pass
1656#ifdef PRINT_RATIO_PROGRESS
1657    iPass++;
1658    printf("Iteration %d pass %d dualOut %g thru %g increase %g numberRemain %d\n",
1659           numberIterations_,iPass,fabs(dualOut_),totalThru+thruThis,increaseInObjective+increaseInThis,
1660           numberRemaining);
1661#endif
1662    if (totalThru + thruThis >= fabs(dualOut_) ||
1663        (increaseInObjective + increaseInThis < 0.0&&bestSequence>=0) || !numberRemaining) {
1664      // We should be pivoting in this batch
1665      // so compress down to this lot
1666      numberRemaining=0;
1667      candidateList.clearAndReset();
1668      for (int iBlock=0;iBlock<numberBlocks;iBlock++) {
1669        int offset=flipList.startPartition(iBlock);
1670        // set so can compact
1671        int numberKeep=result2[iBlock].numberLastSwapped;
1672        flipList.setNumElementsPartition(iBlock,numberKeep);
1673        for (int i = offset+numberKeep; i <offset+result2[iBlock].numberSwapped; i++) {
1674          array[numberRemaining] = flippedElement[i];
1675          flippedElement[i]=0;
1676          indices[numberRemaining++] = flippedIndex[i];
1677        }
1678      }
1679      // make sure candidateList will be cleared correctly
1680      candidateList.setNumElements(numberRemaining);
1681      candidateList.setPackedMode(true);
1682      flipList.compact();
1683      numberSwapped=numberLastSwapped;
1684      useThru=totalThru;
1685      assert (totalThru<=fabs(dualOut_));
1686      int iTry;
1687      // first get ratio with tolerance
1688      for (iTry = 0; iTry < MAXTRY; iTry++) {
1689       
1690        upperTheta = 1.0e50;
1691        numberInList=numberRemaining;
1692        numberRemaining = 0;
1693        increaseInThis = 0.0; //objective increase in this loop
1694        thruThis = 0.0;
1695        for (int i = 0; i < numberInList; i++) {
1696          int iSequence = indices[i];
1697          assert (getInternalStatus(iSequence) != isFree
1698                  && getInternalStatus(iSequence) != superBasic);
1699          int iStatus = internalStatus[iSequence] & 3;
1700          // treat as if at lower bound
1701          double mult = multiplier[iStatus];
1702          double alpha = array[i];
1703          double oldValue = abcDj[iSequence]*mult;
1704          double value = oldValue - upperTheta * alpha;
1705          assert (alpha>0.0);
1706          //double perturbedTolerance = abcPerturbation[iSequence];
1707          if (value < -currentDualTolerance_) {
1708            if (alpha >= currentAcceptablePivot_) {
1709              upperTheta = (oldValue + currentDualTolerance_) / alpha;
1710            }
1711          }
1712        }
1713        bestPivot = currentAcceptablePivot_;
1714        sequenceIn = -1;
1715        double largestPivot = currentAcceptablePivot_;
1716        // Sum of bad small pivots
1717#ifdef MORE_CAREFUL
1718        double sumBadPivots = 0.0;
1719        badSumPivots = false;
1720#endif
1721        // Make sure upperTheta will work (-O2 and above gives problems)
1722        upperTheta *= 1.0000000001;
1723        for (int i = 0; i < numberInList; i++) {
1724          int iSequence = indices[i];
1725          int iStatus = internalStatus[iSequence] & 3;
1726          // treat as if at lower bound
1727          double mult = multiplier[iStatus];
1728          double alpha = array[i];
1729          double oldValue = abcDj[iSequence]*mult;
1730          double value = oldValue - upperTheta * alpha;
1731          assert (alpha>0.0);
1732          //double perturbedTolerance = abcPerturbation[iSequence];
1733          double badDj = 0.0;
1734          if (value <= 0.0) {
1735            // add to list of swapped
1736            flippedElement[numberSwapped] = alpha;
1737            flippedIndex[numberSwapped++] = iSequence;
1738            badDj = oldValue - useTolerance2();
1739            // select if largest pivot
1740            bool take = (alpha > bestPivot); 
1741#ifdef MORE_CAREFUL
1742            if (alpha < currentAcceptablePivot_ && upperTheta < 1.0e20) {
1743              if (value < -currentDualTolerance_) {
1744                double gap = abcUpper[iSequence] - abcLower[iSequence];
1745                if (gap < 1.0e20)
1746                  sumBadPivots -= value * gap;
1747                else
1748                  sumBadPivots += 1.0e20;
1749                //printf("bad %d alpha %g dj at lower %g\n",
1750                //     iSequence,alpha,value);
1751              }
1752            }
1753#endif
1754            if (take&&flagged(iSequence)) {
1755              if (bestPivot>currentAcceptablePivot_)
1756                take=false;
1757            }
1758            if (take) {
1759              sequenceIn = iSequence;
1760              thisPivotValue=alpha;
1761              bestPivot =  alpha;
1762              if (flagged(iSequence)) 
1763                bestPivot=currentAcceptablePivot_;
1764              bestSequence = iSequence;
1765              theta = (oldValue+goToTolerance1()) / alpha;
1766              //assert (oldValue==abcDj[iSequence]);
1767              largestPivot = CoinMax(largestPivot, 0.5 * bestPivot);
1768            }
1769            double range = abcUpper[iSequence] - abcLower[iSequence];
1770            thruThis += range * fabs(alpha);
1771            increaseInThis += badDj * range;
1772          } else {
1773            // add to list of remaining
1774            array[numberRemaining] = alpha;
1775            indices[numberRemaining++] = iSequence;
1776          }
1777        }
1778#ifdef MORE_CAREFUL
1779        // If we have done pivots and things look bad set alpha_ 0.0 to force factorization
1780        if (sumBadPivots > 1.0e4) {
1781#if ABC_NORMAL_DEBUG>0
1782          if (handler_->logLevel() > 1)
1783            printf("maybe forcing re-factorization - sum %g  %d pivots\n", sumBadPivots,
1784                   abcFactorization_->pivots());
1785#endif
1786          if(abcFactorization_->pivots() > 3) {
1787            badSumPivots = true;
1788#ifdef PRINT_RATIO_PROGRESS
1789            printf("maybe forcing re-factorization - sum %g  %d pivots\n", sumBadPivots,
1790                   abcFactorization_->pivots());
1791#endif
1792            break;
1793          }
1794        }
1795#endif
1796        // If we stop now this will be increase in objective (I think)
1797        double increase = (fabs(dualOut_) - totalThru) * theta;
1798        increase += increaseInObjective;
1799        if (theta < 0.0)
1800          thruThis += fabs(dualOut_); // force using this one
1801        if (increaseInObjective < 0.0 && increase < 0.0 && lastSequence >= 0) {
1802          // back
1803          // We may need to be more careful - we could do by
1804          // switch so we always do fine grained?
1805          bestPivot = 0.0;
1806          bestSequence=-1;
1807        } else {
1808          // add in
1809          totalThru += thruThis;
1810          increaseInObjective += increaseInThis;
1811        }
1812        if (bestPivot < 0.1 * bestEverPivot &&
1813            bestEverPivot > 1.0e-6 &&
1814            (bestPivot < 1.0e-3 || totalThru * 2.0 > fabs(dualOut_))) {
1815          // back to previous one
1816          sequenceIn = lastSequence;
1817          thisPivotValue=lastPivotValue;
1818#ifdef PRINT_RATIO_PROGRESS
1819          printf("Try %d back to previous bestPivot %g bestEver %g totalThru %g\n",
1820                 iTry,bestPivot,bestEverPivot,totalThru);
1821#endif
1822          break;
1823        } else if (sequenceIn == -1 && upperTheta > largeValue_) {
1824          if (lastPivotValue > currentAcceptablePivot_) {
1825            // back to previous one
1826            sequenceIn = lastSequence;
1827            thisPivotValue=lastPivotValue;
1828#ifdef PRINT_RATIO_PROGRESS
1829            printf("Try %d no pivot this time large theta %g lastPivot %g\n",
1830                   iTry,upperTheta,lastPivotValue);
1831#endif
1832          } else {
1833            // can only get here if all pivots too small
1834#ifdef PRINT_RATIO_PROGRESS
1835            printf("XXBAD Try %d no pivot this time large theta %g lastPivot %g\n",
1836                   iTry,upperTheta,lastPivotValue);
1837#endif
1838          }
1839          break;
1840        } else if (totalThru >= fabs(dualOut_)) {
1841#ifdef PRINT_RATIO_PROGRESS
1842          printf("Try %d upperTheta %g totalThru %g - can modify costs\n",
1843                 iTry,upperTheta,lastPivotValue);
1844#endif
1845          modifyCosts = true; // fine grain - we can modify costs
1846          break; // no point trying another loop
1847        } else {
1848          lastSequence = sequenceIn;
1849          lastPivotValue=thisPivotValue;
1850          if (bestPivot > bestEverPivot) {
1851            bestEverPivot = bestPivot;
1852            //bestEverSequence=bestSequence;
1853          }
1854          modifyCosts = true; // fine grain - we can modify costs
1855        }
1856#ifdef PRINT_RATIO_PROGRESS
1857        printf("Onto next pass totalthru %g bestPivot %g\n",
1858               totalThru,lastPivotValue);
1859#endif
1860        numberLastSwapped=numberSwapped;
1861      }
1862      break;
1863    } else {
1864      // skip this lot
1865      if (bestPivot > 1.0e-3 || bestPivot > bestEverPivot) {
1866#ifdef PRINT_RATIO_PROGRESS
1867        printf("Continuing bestPivot %g bestEver %g %d swapped - new tentative %g\n",
1868               bestPivot,bestEverPivot,numberSwapped,2*upperTheta);
1869#endif
1870        bestEverPivot = bestPivot;
1871        //bestEverSequence=bestSequence;
1872        lastPivotValue=bestPivot;
1873        lastSequence = bestSequence;
1874        for (int iBlock=0;iBlock<numberBlocks;iBlock++) 
1875          result2[iBlock].numberLastSwapped=result2[iBlock].numberSwapped;
1876        numberLastSwapped=numberSwapped;
1877      } else {
1878#ifdef PRINT_RATIO_PROGRESS
1879        printf("keep old swapped bestPivot %g bestEver %g %d swapped - new tentative %g\n",
1880               bestPivot,bestEverPivot,numberSwapped,2*upperTheta);
1881#endif
1882        // keep old swapped
1883        numberRemaining=0;
1884        for (int iBlock=0;iBlock<numberBlocks;iBlock++) {
1885          int offset=flipList.startPartition(iBlock);
1886          int numberRemaining2=offset+result2[iBlock].numberRemaining;
1887          int numberLastSwapped=result2[iBlock].numberLastSwapped;
1888          for (int i = offset+numberLastSwapped; 
1889               i <offset+result2[iBlock].numberSwapped; i++) {
1890            array[numberRemaining2] = flippedElement[i];
1891            flippedElement[i]=0;
1892            indices[numberRemaining2++] = flippedIndex[i];
1893          }
1894          result2[iBlock].numberRemaining=numberRemaining2-offset;
1895          numberRemaining += result2[iBlock].numberRemaining;
1896          result2[iBlock].numberSwapped=numberLastSwapped;
1897        }
1898        numberSwapped=numberLastSwapped;
1899      }
1900      increaseInObjective += increaseInThis;
1901      tentativeTheta = 2.0 * upperTheta;
1902      totalThru += thruThis;
1903    }
1904  }
1905  if (flipList.getNumPartitions()) {
1906    // need to compress
1907    numberRemaining=0;
1908    candidateList.clearAndReset();
1909    for (int iBlock=0;iBlock<numberBlocks;iBlock++) {
1910      int offset=flipList.startPartition(iBlock);
1911      // set so can compact
1912      int numberKeep=result2[iBlock].numberLastSwapped;
1913      flipList.setNumElementsPartition(iBlock,numberKeep);
1914      for (int i = offset+numberKeep; i <offset+result2[iBlock].numberSwapped; i++) {
1915        array[numberRemaining] = flippedElement[i];
1916        flippedElement[i]=0;
1917        indices[numberRemaining++] = flippedIndex[i];
1918      }
1919    }
1920    // make sure candidateList will be cleared correctly
1921    candidateList.setNumElements(numberRemaining);
1922    candidateList.setPackedMode(true);
1923    flipList.compact();
1924    numberSwapped=numberLastSwapped;
1925  }
1926  result.theta=theta; //?
1927  result.totalThru=totalThru;
1928  result.useThru=useThru;
1929  result.bestEverPivot=bestEverPivot; //?
1930  result.increaseInObjective=increaseInObjective; //?
1931  result.tentativeTheta=tentativeTheta; //?
1932  result.lastPivotValue=lastPivotValue;
1933  result.thisPivotValue=thisPivotValue;
1934  result.lastSequence=lastSequence;
1935  result.sequence=sequenceIn;
1936  //result.block;
1937  //result.pass;
1938  //result.phase;
1939  result.numberSwapped=numberSwapped;
1940  result.numberLastSwapped=numberLastSwapped;
1941  result.modifyCosts=modifyCosts;
1942}
1943/*
1944  Chooses incoming
1945  Puts flipped ones in list
1946  If necessary will modify costs
1947*/
1948void
1949AbcSimplexDual::dualColumn2()
1950{
1951  CoinPartitionedVector & candidateList = usefulArray_[arrayForDualColumn_]; 
1952  CoinPartitionedVector & flipList = usefulArray_[arrayForFlipBounds_];
1953  //usefulArray_[arrayForTableauRow_].compact();
1954  usefulArray_[arrayForTableauRow_].computeNumberElements();
1955  if (candidateList.getNumElements()<100)
1956    candidateList.compact();
1957  //usefulArray_[arrayForTableauRow_].compact();
1958  // But factorizations complain if <1.0e-8
1959  //acceptablePivot=CoinMax(acceptablePivot,1.0e-8);
1960  // do ratio test for normal iteration
1961  currentAcceptablePivot_ = 1.0e-1 * acceptablePivot_;
1962  if (numberIterations_ > 100)
1963    currentAcceptablePivot_ = acceptablePivot_;
1964  if (abcFactorization_->pivots() > 10 ||
1965      (abcFactorization_->pivots() && sumDualInfeasibilities_))
1966    currentAcceptablePivot_ = 1.0e+3 * acceptablePivot_; // if we have iterated be more strict
1967  else if (abcFactorization_->pivots() > 5)
1968    currentAcceptablePivot_ = 1.0e+2 * acceptablePivot_; // if we have iterated be slightly more strict
1969  else if (abcFactorization_->pivots())
1970    currentAcceptablePivot_ = acceptablePivot_; // relax
1971  // If very large infeasibility use larger acceptable pivot
1972  if (abcFactorization_->pivots()) {
1973    double value=CoinMax(currentAcceptablePivot_,1.0e-12*fabs(dualOut_));
1974    currentAcceptablePivot_=CoinMin(1.0e2*currentAcceptablePivot_,value);
1975  }
1976  // return at once if no free but there were free
1977  if (lastFirstFree_>=0&&sequenceIn_<0) {
1978    candidateList.clearAndReset();
1979    upperTheta_=COIN_DBL_MAX;
1980  }
1981  if (upperTheta_==COIN_DBL_MAX) {
1982    usefulArray_[arrayForTableauRow_].compact();
1983    assert (!candidateList.getNumElements());
1984    return;
1985  }
1986  int numberSwapped = 0;
1987  //int numberRemaining = candidateList.getNumElements();
1988 
1989  double useThru = 0.0; // for when variables flip
1990 
1991  // If we think we need to modify costs (not if something from broad sweep)
1992  bool modifyCosts = false;
1993  // pivot elements
1994  //double *  COIN_RESTRICT array=candidateList.denseVector();
1995  // indices
1996  const double *  COIN_RESTRICT abcLower = abcLower_;
1997  const double *  COIN_RESTRICT abcUpper = abcUpper_;
1998  const double *  COIN_RESTRICT abcSolution = abcSolution_;
1999  const double *  COIN_RESTRICT abcDj = abcDj_;
2000  //const double *  COIN_RESTRICT abcPerturbation = abcPerturbation_;
2001  const unsigned char *  COIN_RESTRICT internalStatus = internalStatus_;
2002  // We can also see if infeasible or pivoting on free
2003  const double multiplier[] = { 1.0, -1.0};
2004
2005  // For list of flipped
2006  int numberLastSwapped=0;
2007  int *  COIN_RESTRICT flippedIndex = flipList.getIndices();
2008  double *  COIN_RESTRICT flippedElement = flipList.denseVector();
2009  // If sum of bad small pivots too much
2010#ifdef MORE_CAREFUL
2011  bool badSumPivots = false;
2012#endif
2013  if (sequenceIn_<0) {
2014    int lastSequence = -1;
2015    double lastPivotValue = 0.0;
2016    double thisPivotValue=0.0;
2017    dualColumnResult result;
2018    dualColumn2Most(result);
2019    sequenceIn_=result.sequence;
2020    lastSequence=result.lastSequence;
2021    thisPivotValue=result.thisPivotValue;
2022    lastPivotValue=result.lastPivotValue;
2023    numberSwapped=result.numberSwapped;
2024    numberLastSwapped=result.numberLastSwapped;
2025    modifyCosts=result.modifyCosts;
2026    useThru=result.useThru;
2027#ifdef PRINT_RATIO_PROGRESS
2028    printf("Iteration %d - out of loop - sequenceIn %d pivotValue %g\n",
2029           numberIterations_,sequenceIn_,thisPivotValue);
2030#endif
2031    // can get here without sequenceIn_ set but with lastSequence
2032    if (sequenceIn_ < 0 && lastSequence >= 0) {
2033      // back to previous one
2034      sequenceIn_ = lastSequence;
2035      thisPivotValue=lastPivotValue;
2036#ifdef PRINT_RATIO_PROGRESS
2037      printf("BUT back one - out of loop - sequenceIn %d pivotValue %g\n",
2038           sequenceIn_,thisPivotValue);
2039#endif
2040    }
2041   
2042    // Movement should be minimum for anti-degeneracy - unless
2043    // fixed variable out
2044    // no need if just one candidate
2045    assert( numberSwapped-numberLastSwapped>=0);
2046    double minimumTheta;
2047    if (upperOut_ > lowerOut_&&numberSwapped-numberLastSwapped>1)
2048#ifndef TRY_ABC_GUS
2049      minimumTheta = MINIMUMTHETA;
2050#else
2051      minimumTheta = minimumThetaMovement_;
2052#endif
2053    else
2054      minimumTheta = 0.0;
2055    if (sequenceIn_ >= 0) {
2056      double oldValue;
2057      alpha_ = thisPivotValue; 
2058      // correct sign
2059      int iStatus = internalStatus[sequenceIn_] & 3;
2060      alpha_ *= multiplier[iStatus];
2061      oldValue = abcDj[sequenceIn_];
2062      theta_ = CoinMax(oldValue / alpha_, 0.0);
2063      if (theta_ < minimumTheta && fabs(alpha_) < 1.0e5 && 1) {
2064        // can't pivot to zero
2065        theta_ = minimumTheta;
2066      }
2067      // may need to adjust costs so all dual feasible AND pivoted is exactly 0
2068      if (modifyCosts
2069#ifdef MORE_CAREFUL
2070          &&!badSumPivots
2071#endif
2072          ) {
2073        for (int i = numberLastSwapped; i < numberSwapped; i++) {
2074          int iSequence = flippedIndex[i];
2075          int iStatus = internalStatus[iSequence] & 3;
2076          // treat as if at lower bound
2077          double mult = multiplier[iStatus];
2078          double alpha = flippedElement[i];
2079          flippedElement[i]=0.0;
2080          if (iSequence==sequenceIn_)
2081            continue;
2082          double oldValue = abcDj[iSequence]*mult;
2083          double value = oldValue - theta_ * alpha;
2084          //double perturbedTolerance = abcPerturbation[iSequence];
2085          if (-value > currentDualTolerance_) {
2086            //thisIncrease = true;
2087#if MODIFYCOST
2088            if (numberIterations_>=normalDualColumnIteration_) {
2089              // modify cost to hit new tolerance
2090              double modification = alpha * theta_ - oldValue - currentDualTolerance_;
2091              if ((specialOptions_&(2048 + 4096)) != 0) {
2092                if ((specialOptions_ & 2048) != 0) {
2093                  if (fabs(modification) < 1.0e-10)
2094                    modification = 0.0;
2095                } else {
2096                  if (fabs(modification) < 1.0e-12)
2097                    modification = 0.0;
2098                }
2099              }
2100              abcDj_[iSequence] += mult*modification;
2101              abcCost_[iSequence] +=  mult*modification;
2102              if (modification)
2103                numberChanged_ ++; // Say changed costs
2104            }
2105#endif
2106          }
2107        }
2108        numberSwapped=numberLastSwapped;
2109      }
2110    }
2111  }
2112 
2113#ifdef MORE_CAREFUL
2114  // If we have done pivots and things look bad set alpha_ 0.0 to force factorization
2115  if ((badSumPivots ||
2116       fabs(theta_ * badFree) > 10.0 * currentDualTolerance_) && abcFactorization_->pivots()) {
2117#if ABC_NORMAL_DEBUG>0
2118    if (handler_->logLevel() > 1)
2119      printf("forcing re-factorization\n");
2120#endif
2121    sequenceIn_ = -1;
2122  }
2123#endif
2124  if (sequenceIn_ >= 0) {
2125    lowerIn_ = abcLower[sequenceIn_];
2126    upperIn_ = abcUpper[sequenceIn_];
2127    valueIn_ = abcSolution[sequenceIn_];
2128    dualIn_ = abcDj_[sequenceIn_];
2129   
2130    if (numberTimesOptimal_) {
2131      // can we adjust cost back closer to original
2132      //*** add coding
2133    }
2134    if (numberIterations_>=normalDualColumnIteration_) {
2135#if MODIFYCOST>1
2136      // modify cost to hit zero exactly
2137      // so (dualIn_+modification)==theta_*alpha_
2138      //double perturbedTolerance = abcPerturbation[sequenceIn_];
2139      double modification = theta_ * alpha_ - (dualIn_+goToTolerance2());
2140      // But should not move objective too much ??
2141#ifdef DONT_MOVE_OBJECTIVE
2142      double moveObjective = fabs(modification * abcSolution[sequenceIn_]);
2143      double smallMove = CoinMax(fabs(rawObjectiveValue_), 1.0e-3);
2144      if (moveObjective > smallMove) {
2145#if ABC_NORMAL_DEBUG>0
2146        if (handler_->logLevel() > 1)
2147          printf("would move objective by %g - original mod %g sol value %g\n", moveObjective,
2148                 modification, abcSolution[sequenceIn_]);
2149#endif
2150        modification *= smallMove / moveObjective;
2151      }
2152#endif
2153#ifdef MORE_CAREFUL
2154      if (badSumPivots)
2155        modification = 0.0;
2156#endif
2157      if (atFakeBound(sequenceIn_)&&fabs(modification)<currentDualTolerance_)
2158        modification=0.0;
2159      if ((specialOptions_&(2048 + 4096)) != 0) {
2160        if ((specialOptions_ & 16384) != 0) {
2161          // in fast dual
2162          if (fabs(modification) < 1.0e-7)
2163            modification = 0.0;
2164        } else if ((specialOptions_ & 2048) != 0) {
2165          if (fabs(modification) < 1.0e-10)
2166            modification = 0.0;
2167        } else {
2168          if (fabs(modification) < 1.0e-12)
2169            modification = 0.0;
2170        }
2171      }
2172     
2173      dualIn_ += modification;
2174      abcDj_[sequenceIn_] = dualIn_;
2175      abcCost_[sequenceIn_] += modification;
2176      if (modification)
2177        numberChanged_ ++; // Say changed costs
2178      //int costOffset = numberRows_+numberColumns_;
2179      //abcCost[sequenceIn_+costOffset] += modification; // save change
2180      //assert (fabs(modification)<1.0e-6);
2181#if ABC_NORMAL_DEBUG>3
2182      if ((handler_->logLevel() & 32) && fabs(modification) > 1.0e-15)
2183        printf("exact %d new cost %g, change %g\n", sequenceIn_,
2184               abcCost_[sequenceIn_], modification);
2185#endif
2186#endif
2187    }
2188   
2189    if (alpha_ < 0.0) {
2190      // as if from upper bound
2191      directionIn_ = -1;
2192      upperIn_ = valueIn_;
2193    } else {
2194      // as if from lower bound
2195      directionIn_ = 1;
2196      lowerIn_ = valueIn_;
2197    }
2198  } else {
2199    // no pivot
2200    alpha_ = 0.0;
2201  }
2202  // clear array
2203  //printf("dualout %g theta %g alpha %g dj %g thru %g swapped %d product %g guess %g useThru %g\n",
2204  //     dualOut_,theta_,alpha_,abcDj_[sequenceIn_],useThru,numberSwapped,
2205  //     fabs(dualOut_)*theta_,objectiveValue()+(fabs(dualOut_)-useThru)*theta_,useThru);
2206  objectiveChange_ = (fabs(dualOut_)-useThru)*theta_;
2207  rawObjectiveValue_+= objectiveChange_;
2208  //CoinZeroN(array, maximumNumberInList);
2209  //candidateList.setNumElements(0);
2210  candidateList.clearAndReset();
2211  flipList.setNumElements(numberSwapped);
2212  if (numberSwapped)
2213    flipList.setPackedMode(true);
2214#ifdef PRINT_RATIO_PROGRESS
2215  printf("%d flipped\n",numberSwapped);
2216#endif
2217}
2218/* Checks if tentative optimal actually means unbounded
2219   Returns -3 if not, 2 if is unbounded */
2220int
2221AbcSimplexDual::checkUnbounded(CoinIndexedVector & ray,
2222                               double changeCost)
2223{
2224  int status = 2; // say unbounded
2225  abcFactorization_->updateColumn(ray);
2226  // get reduced cost
2227  int number = ray.getNumElements();
2228  int *  COIN_RESTRICT index = ray.getIndices();
2229  double *  COIN_RESTRICT array = ray.denseVector();
2230  const int * COIN_RESTRICT abcPivotVariable = abcPivotVariable_;
2231  for (int i = 0; i < number; i++) {
2232    int iRow = index[i];
2233    int iPivot = abcPivotVariable[iRow];
2234    changeCost -= cost(iPivot) * array[iRow];
2235  }
2236  double way;
2237  if (changeCost > 0.0) {
2238    //try going down
2239    way = 1.0;
2240  } else if (changeCost < 0.0) {
2241    //try going up
2242    way = -1.0;
2243  } else {
2244#if ABC_NORMAL_DEBUG>3
2245    printf("can't decide on up or down\n");
2246#endif
2247    way = 0.0;
2248    status = -3;
2249  }
2250  double movement = 1.0e10 * way; // some largish number
2251  double zeroTolerance = 1.0e-14 * dualBound_;
2252  for (int i = 0; i < number; i++) {
2253    int iRow = index[i];
2254    int iPivot = abcPivotVariable[iRow];
2255    double arrayValue = array[iRow];
2256    if (fabs(arrayValue) < zeroTolerance)
2257      arrayValue = 0.0;
2258    double newValue = solution(iPivot) + movement * arrayValue;
2259    if (newValue > upper(iPivot) + primalTolerance_ ||
2260        newValue < lower(iPivot) - primalTolerance_)
2261      status = -3; // not unbounded
2262  }
2263  if (status == 2) {
2264    // create ray
2265    delete [] ray_;
2266    ray_ = new double [numberColumns_];
2267    CoinZeroN(ray_, numberColumns_);
2268    for (int i = 0; i < number; i++) {
2269      int iRow = index[i];
2270      int iPivot = abcPivotVariable[iRow];
2271      double arrayValue = array[iRow];
2272      if (iPivot < numberColumns_ && fabs(arrayValue) >= zeroTolerance)
2273        ray_[iPivot] = way * array[iRow];
2274    }
2275  }
2276  ray.clear();
2277  return status;
2278}
2279// Saves good status etc
2280void 
2281AbcSimplex::saveGoodStatus()
2282{
2283  // save arrays
2284  CoinMemcpyN(internalStatus_, numberTotal_, internalStatusSaved_);
2285  CoinMemcpyN(abcSolution_,numberTotal_,solutionSaved_);
2286}
2287// Restores previous good status and says trouble
2288void 
2289AbcSimplex::restoreGoodStatus(int type)
2290{
2291#if PARTITION_ROW_COPY==1
2292  stateOfProblem_ |= NEED_BASIS_SORT;
2293#endif
2294  if (type) {
2295    if (alphaAccuracy_ != -1.0)
2296      alphaAccuracy_ = -2.0;
2297    // trouble - restore solution
2298    CoinMemcpyN(internalStatusSaved_, numberTotal_, internalStatus_);
2299    CoinMemcpyN(solutionSaved_,numberTotal_, abcSolution_);
2300    int *  COIN_RESTRICT abcPivotVariable = abcPivotVariable_;
2301    // redo pivotVariable
2302    int numberBasic=0;
2303    for (int i=0;i<numberTotal_;i++) {
2304      if (getInternalStatus(i)==basic)
2305        abcPivotVariable[numberBasic++]=i;
2306    }
2307    assert (numberBasic==numberRows_);
2308    forceFactorization_ = 1; // a bit drastic but ..
2309    changeMade_++; // say something changed
2310  }
2311}
2312static int computePrimalsAndCheck(AbcSimplexDual * model,const int whichArray[2])
2313{
2314  CoinIndexedVector * array1 = model->usefulArray(whichArray[0]);
2315  CoinIndexedVector * array2 = model->usefulArray(whichArray[1]);
2316  int numberRefinements=model->computePrimals(array1,array2);
2317  model->checkPrimalSolution(true);
2318  return numberRefinements;
2319}
2320static int computeDualsAndCheck(AbcSimplexDual * model,const int whichArray[2])
2321{
2322  CoinIndexedVector * array1 = model->usefulArray(whichArray[0]);
2323  CoinIndexedVector * array2 = model->usefulArray(whichArray[1]);
2324  int numberRefinements=model->computeDuals(NULL,array1,array2);
2325  model->checkDualSolutionPlusFake();
2326  return numberRefinements;
2327}
2328// Computes solutions - 1 do duals, 2 do primals, 3 both
2329int
2330AbcSimplex::gutsOfSolution(int type)
2331{
2332  AbcSimplexDual * dual = reinterpret_cast<AbcSimplexDual *>(this);
2333  // do work and check
2334  int numberRefinements=0;
2335#if ABC_PARALLEL
2336  if (type!=3) {
2337#endif
2338    if ((type&2)!=0) {
2339      //work space
2340      int whichArray[2];
2341      for (int i=0;i<2;i++)
2342        whichArray[i]=getAvailableArray();
2343      numberRefinements=computePrimalsAndCheck(dual,whichArray);
2344      for (int i=0;i<2;i++)
2345        setAvailableArray(whichArray[i]);
2346    }
2347    if ((type&1)!=0
2348#if CAN_HAVE_ZERO_OBJ>1
2349        &&(specialOptions_&2097152)==0
2350#endif
2351        ) {
2352      //work space
2353      int whichArray[2];
2354      for (int i=0;i<2;i++)
2355        whichArray[i]=getAvailableArray();
2356      numberRefinements+=computeDualsAndCheck(dual,whichArray);
2357      for (int i=0;i<2;i++)
2358        setAvailableArray(whichArray[i]);
2359    }
2360#if ABC_PARALLEL
2361  } else {
2362#ifndef NDEBUG
2363    abcFactorization_->checkMarkArrays();
2364#endif
2365    // redo to do in parallel
2366      //work space
2367      int whichArray[5];
2368      for (int i=1;i<5;i++)
2369        whichArray[i]=getAvailableArray();
2370      if (parallelMode_==0) {
2371        numberRefinements+=computePrimalsAndCheck(dual,whichArray+3);
2372        numberRefinements+=computeDualsAndCheck(dual,whichArray+1);
2373      } else {
2374#if ABC_PARALLEL==1
2375        threadInfo_[0].stuff[0]=whichArray[1];
2376        threadInfo_[0].stuff[1]=whichArray[2];
2377        stopStart_=1+32*1;
2378        startParallelStuff(1);
2379#else
2380        whichArray[0]=1;
2381        CoinAbcThreadInfo info;
2382        info.status=1;
2383        info.stuff[0]=whichArray[1];
2384        info.stuff[1]=whichArray[2];
2385        int n=cilk_spawn computeDualsAndCheck(dual,whichArray+1);
2386#endif
2387        numberRefinements=computePrimalsAndCheck(dual,whichArray+3);
2388#if ABC_PARALLEL==1
2389        numberRefinements+=stopParallelStuff(1);
2390#else
2391        cilk_sync;
2392        numberRefinements+=n;
2393#endif
2394      }
2395      for (int i=1;i<5;i++)
2396        setAvailableArray(whichArray[i]);
2397  }
2398#endif
2399  rawObjectiveValue_ +=sumNonBasicCosts_;
2400  setClpSimplexObjectiveValue();
2401  if (handler_->logLevel() > 3 || (largestPrimalError_ > 1.0e-2 ||
2402                                   largestDualError_ > 1.0e-2))
2403    handler_->message(CLP_SIMPLEX_ACCURACY, messages_)
2404      << largestPrimalError_
2405      << largestDualError_
2406      << CoinMessageEol;
2407  if ((largestPrimalError_ > 1.0e1||largestDualError_>1.0e1) 
2408      && numberRows_ > 100 && numberIterations_) {
2409#if ABC_NORMAL_DEBUG>0
2410    printf("Large errors - primal %g dual %g\n",
2411           largestPrimalError_,largestDualError_);
2412#endif
2413    // Change factorization tolerance
2414    //if (abcFactorization_->zeroTolerance() > 1.0e-18)
2415    //abcFactorization_->zeroTolerance(1.0e-18);
2416  }
2417  return numberRefinements;
2418}
2419// Make non free variables dual feasible by moving to a bound
2420int 
2421AbcSimplexDual::makeNonFreeVariablesDualFeasible(bool changeCosts)
2422{
2423  const double multiplier[] = { 1.0, -1.0}; // sign?????????????
2424  double dualT = - 100.0*currentDualTolerance_;
2425  int numberMoved=0;
2426  int numberChanged=0;
2427  double *  COIN_RESTRICT solution = abcSolution_;
2428  double *  COIN_RESTRICT lower = abcLower_;
2429  double *  COIN_RESTRICT upper = abcUpper_;
2430  double saveNonBasicCosts =sumNonBasicCosts_;
2431  double largestCost=dualTolerance_;
2432  for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
2433    unsigned char iStatus = internalStatus_[iSequence] & 7;
2434    if (iStatus<4) {
2435      double mult = multiplier[iStatus];
2436      double dj = abcDj_[iSequence] * mult;
2437      if (dj < dualT) {
2438        largestCost=CoinMax(fabs(abcPerturbation_[iSequence]),largestCost);
2439      }
2440    } else if (iStatus!=7) {
2441      largestCost=CoinMax(fabs(abcPerturbation_[iSequence]),largestCost);
2442    }
2443  }
2444  for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
2445    unsigned char iStatus = internalStatus_[iSequence] & 7;
2446    if (iStatus<4) {
2447      double mult = multiplier[iStatus];
2448      double dj = abcDj_[iSequence] * mult;
2449      // ? should use perturbationArray
2450      if (dj < dualT&&changeCosts) {
2451        double gap=upper[iSequence]-lower[iSequence];
2452        if (gap>1.1*currentDualBound_) {
2453          assert(getFakeBound(iSequence)==AbcSimplexDual::noFake);
2454          if (!iStatus) {
2455            upper[iSequence]=lower[iSequence]+currentDualBound_;
2456            setFakeBound(iSequence, AbcSimplexDual::upperFake);
2457          } else {
2458            lower[iSequence]=upper[iSequence]-currentDualBound_;
2459            setFakeBound(iSequence, AbcSimplexDual::lowerFake);
2460          }
2461        }
2462        double costDifference=abcPerturbation_[iSequence]-abcCost_[iSequence];
2463        if (costDifference*mult>0.0) {
2464#if ABC_NORMAL_DEBUG>0
2465          //printf("can improve situation on %d by %g\n",iSequence,costDifference*mult);
2466#endif
2467          dj+=costDifference*mult;
2468          abcDj_[iSequence]+=costDifference;
2469          abcCost_[iSequence]=abcPerturbation_[iSequence];
2470          costDifference=0.0;
2471        }
2472        double changeInObjective=-gap*dj;
2473        // probably need to be more intelligent
2474        if ((changeInObjective>1.0e3||gap>CoinMax(0.5*dualBound_,1.0e3))&&dj>-1.0e-1) {
2475          // change cost
2476          costDifference=(-dj-0.5*dualTolerance_)*mult;
2477          if (fabs(costDifference)<100000.0*dualTolerance_) {
2478            dj+=costDifference*mult;
2479            abcDj_[iSequence]+=costDifference;
2480            abcCost_[iSequence]+=costDifference;
2481            numberChanged++;
2482          }
2483        }
2484      }
2485      if (dj < dualT&&!changeCosts) {
2486        numberMoved++;
2487        if (iStatus==0) {
2488          // at lower bound
2489          // to upper bound
2490          setInternalStatus(iSequence, atUpperBound);
2491          solution[iSequence] = upper[iSequence];
2492          sumNonBasicCosts_ += abcCost_[iSequence]*(upper[iSequence]-lower[iSequence]);
2493        } else {
2494          // at upper bound
2495          // to lower bound
2496          setInternalStatus(iSequence, atLowerBound);
2497          solution[iSequence] = lower[iSequence];
2498          sumNonBasicCosts_ -= abcCost_[iSequence]*(upper[iSequence]-lower[iSequence]);
2499        }
2500      }
2501    }
2502  }
2503#if ABC_NORMAL_DEBUG>0
2504  if (numberMoved+numberChanged)
2505  printf("makeDualfeasible %d moved, %d had costs changed - sum dual %g\n",
2506         numberMoved,numberChanged,sumDualInfeasibilities_);
2507#endif
2508  rawObjectiveValue_ +=(sumNonBasicCosts_-saveNonBasicCosts);
2509  setClpSimplexObjectiveValue();
2510  return numberMoved;
2511}
2512int 
2513AbcSimplexDual::whatNext()
2514{
2515  int returnCode=0;
2516  assert (!stateOfIteration_);
2517  // see if update stable
2518#if ABC_NORMAL_DEBUG>3
2519  if ((handler_->logLevel() & 32))
2520    printf("btran alpha %g, ftran alpha %g\n", btranAlpha_, alpha_);
2521#endif
2522  int numberPivots=abcFactorization_->pivots();
2523  //double checkValue = 1.0e-7; // numberPivots<5 ? 1.0e-7 : 1.0e-6;
2524  double checkValue = numberPivots ? 1.0e-7 : 1.0e-5;
2525  // relax if free variable in
2526  if ((internalStatus_[sequenceIn_]&7)>1)
2527    checkValue=1.0e-5;
2528  // if can't trust much and long way from optimal then relax
2529  if (largestPrimalError_ > 10.0)
2530    checkValue = CoinMin(1.0e-4, 1.0e-8 * largestPrimalError_);
2531  if (fabs(btranAlpha_) < 1.0e-12 || fabs(alpha_) < 1.0e-12 ||
2532      fabs(btranAlpha_ - alpha_) > checkValue*(1.0 + fabs(alpha_))) {
2533    handler_->message(CLP_DUAL_CHECK, messages_)
2534      << btranAlpha_
2535      << alpha_
2536      << CoinMessageEol;
2537    // see with more relaxed criterion
2538    double test;
2539    if (fabs(btranAlpha_) < 1.0e-8 || fabs(alpha_) < 1.0e-8)
2540      test = 1.0e-1 * fabs(alpha_);
2541    else
2542      test = 10.0*checkValue;//1.0e-4 * (1.0 + fabs(alpha_));
2543    bool needFlag = (fabs(btranAlpha_) < 1.0e-12 || fabs(alpha_) < 1.0e-12 ||
2544                     fabs(btranAlpha_ - alpha_) > test);
2545#if ABC_NORMAL_DEBUG>0
2546    if (fabs(btranAlpha_ + alpha_) < 1.0e-5*(1.0 + fabs(alpha_))) {
2547      printf("alpha diff (%g,%g) %g check %g\n",btranAlpha_,alpha_,fabs(btranAlpha_-alpha_),checkValue*(1.0+fabs(alpha_)));
2548    }
2549#endif
2550    if (numberPivots) {
2551      if (needFlag&&numberPivots<10) {
2552        // need to reject something
2553        assert (sequenceOut_>=0);
2554        char x = isColumn(sequenceOut_) ? 'C' : 'R';
2555        handler_->message(CLP_SIMPLEX_FLAG, messages_)
2556          << x << sequenceWithin(sequenceOut_)
2557          << CoinMessageEol;
2558#if ABC_NORMAL_DEBUG>0
2559        printf("flag %d as alpha  %g %g - after %d pivots\n", sequenceOut_,
2560               btranAlpha_, alpha_,numberPivots);
2561#endif
2562        // Make safer?
2563        abcFactorization_->saferTolerances (1.0, -1.03);
2564        setFlagged(sequenceOut_);
2565        // so can't happen immediately again
2566        sequenceOut_=-1;
2567        lastBadIteration_ = numberIterations_; // say be more cautious
2568      }
2569      // Make safer?
2570      //if (numberPivots<5)
2571      //abcFactorization_->saferTolerances (-0.99, -1.01);
2572      problemStatus_ = -2; // factorize now
2573      returnCode = -2;
2574      stateOfIteration_=2;
2575    } else {
2576      if (needFlag) {
2577        // need to reject something
2578        assert (sequenceOut_>=0);
2579        char x = isColumn(sequenceOut_) ? 'C' : 'R';
2580        handler_->message(CLP_SIMPLEX_FLAG, messages_)
2581          << x << sequenceWithin(sequenceOut_)
2582          << CoinMessageEol;
2583#if ABC_NORMAL_DEBUG>3
2584        printf("flag a %g %g\n", btranAlpha_, alpha_);
2585#endif
2586        // Make safer?
2587        double tolerance=abcFactorization_->pivotTolerance();
2588        abcFactorization_->saferTolerances (1.0, -1.03);
2589        setFlagged(sequenceOut_);
2590        // so can't happen immediately again
2591        sequenceOut_=-1;
2592        //abcProgress_.clearBadTimes();
2593        lastBadIteration_ = numberIterations_; // say be more cautious
2594        if (fabs(alpha_) < 1.0e-10 && fabs(btranAlpha_) < 1.0e-8 && numberIterations_ > 100) {
2595          //printf("I think should declare infeasible\n");
2596          problemStatus_ = 1;
2597          returnCode = 1;
2598          stateOfIteration_=2;
2599        } else {
2600          stateOfIteration_=1;
2601          if (tolerance<abcFactorization_->pivotTolerance()) {
2602            stateOfIteration_=2;
2603            returnCode=-2;
2604          }
2605        }
2606      }
2607    }
2608  }
2609  if (!stateOfIteration_) {
2610    // check update
2611    returnCode=abcFactorization_->checkReplacePart2(pivotRow_,btranAlpha_,alpha_,ftAlpha_);
2612  }
2613  return returnCode;
2614}
2615/* Checks if finished. 
2616   type 0 initial
2617   1 normal
2618   2 initial but we have solved problem before
2619   3 can skip factorization
2620   Updates status */
2621void
2622AbcSimplexDual::statusOfProblemInDual(int type)
2623{
2624  // If lots of iterations then adjust costs if large ones
2625  int numberPivots = abcFactorization_->pivots();
2626#if 0
2627  printf("statusOf %d pivots type %d pstatus %d forceFac %d dont %d\n",
2628         numberPivots,type,problemStatus_,forceFactorization_,dontFactorizePivots_);
2629#endif
2630  int saveType=type;
2631  double lastSumDualInfeasibilities=sumDualInfeasibilities_;
2632  double lastSumPrimalInfeasibilities=sumPrimalInfeasibilities_;
2633  if (type>=4) {
2634    // force factorization
2635    type &= 3;
2636    numberPivots=9999999;
2637  }
2638  //double realDualInfeasibilities = 0.0;
2639  double changeCost;
2640  bool unflagVariables = numberFlagged_>0;
2641  int weightsSaved = 0;
2642  int numberBlackMarks=0;
2643  if (!type) {
2644    // be optimistic
2645    stateOfProblem_ &= ~PESSIMISTIC;
2646    // but use 0.1
2647    double newTolerance = CoinMax(0.1, saveData_.pivotTolerance_);
2648    abcFactorization_->pivotTolerance(newTolerance);
2649  }
2650  if (type != 3&&(numberPivots>dontFactorizePivots_||numberIterations_==baseIteration_)) {
2651    // factorize
2652    // save dual weights
2653    abcDualRowPivot_->saveWeights(this, 1);
2654    weightsSaved = 2;
2655    // is factorization okay?
2656    int factorizationStatus=0;
2657  if ((largestPrimalError_ > 1.0e1||largestDualError_>1.0e1) 
2658      && numberRows_ > 100 && numberIterations_>baseIteration_) {
2659    double newTolerance = CoinMin(1.1* abcFactorization_->pivotTolerance(), 0.99);
2660    abcFactorization_->pivotTolerance(newTolerance);
2661  }
2662#ifdef EARLY_FACTORIZE
2663    if ((numberEarly_&0xffff)<=5||!usefulArray_[ABC_NUMBER_USEFUL-1].getNumElements()) {
2664        //#define KEEP_TRYING_EARLY
2665#ifdef KEEP_TRYING_EARLY
2666      int savedEarly=numberEarly_>>16;
2667      int newEarly=(numberEarly_&0xffff);
2668      newEarly=CoinMin(savedEarly,1+2*newEarly);
2669      numberEarly_=(savedEarly<<16)|newEarly;
2670#endif
2671      factorizationStatus = internalFactorize(type ? 1 : 2);
2672    } else {
2673      CoinIndexedVector & vector = usefulArray_[ABC_NUMBER_USEFUL-1];
2674      bool bad = vector.getNumElements()<0;
2675      const int * indices = vector.getIndices();
2676      double * array = vector.denseVector();
2677      int capacity = vector.capacity();
2678      capacity--;
2679      int numberPivotsStored = indices[capacity];
2680#if ABC_NORMAL_DEBUG>0
2681      printf("early status %s nPivots %d\n",bad ? "bad" : "good",numberPivotsStored);
2682#endif
2683      if (!bad) {
2684        // do remaining pivots
2685        int iterationsDone=abcEarlyFactorization_->pivots();
2686        factorizationStatus = 
2687          abcEarlyFactorization_->replaceColumns(this,vector,
2688                                                 iterationsDone,numberPivotsStored,true);
2689        if (!factorizationStatus) {
2690          // should be able to compare e.g. basics 1.0
2691#define GOOD_EARLY
2692#ifdef GOOD_EARLY
2693          AbcSimplexFactorization * temp = abcFactorization_;
2694          abcFactorization_=abcEarlyFactorization_;
2695          abcEarlyFactorization_=temp;
2696#endif
2697          moveToBasic();
2698#ifndef GOOD_EARLY
2699          int * pivotSave = CoinCopyOfArray(abcPivotVariable_,numberRows_);
2700          factorizationStatus = internalFactorize(type ? 1 : 2);
2701          {
2702            double * array = usefulArray_[3].denseVector();
2703            double * array2 = usefulArray_[4].denseVector();
2704            for (int iPivot = 0; iPivot < numberRows_; iPivot++) {
2705              int iSequence = abcPivotVariable_[iPivot];
2706              unpack(usefulArray_[3], iSequence);
2707              unpack(usefulArray_[4], iSequence);
2708              abcFactorization_->updateColumn(usefulArray_[3]);
2709              abcEarlyFactorization_->updateColumn(usefulArray_[4]);
2710              assert (fabs(array[iPivot] - 1.0) < 1.0e-4);
2711              array[iPivot] = 0.0;
2712              int iPivot2;
2713              for (iPivot2=0;iPivot2<numberRows_;iPivot2++) {
2714                if (pivotSave[iPivot2]==iSequence)
2715                  break;
2716              }
2717              assert (iPivot2<numberRows_);
2718              assert (fabs(array2[iPivot2] - 1.0) < 1.0e-4);
2719              array2[iPivot2] = 0.0;
2720              for (int i = 0; i < numberRows_; i++)
2721                assert (fabs(array[i]) < 1.0e-4);
2722              for (int i = 0; i < numberRows_; i++)
2723                assert (fabs(array2[i]) < 1.0e-4);
2724              usefulArray_[3].clear();
2725              usefulArray_[4].clear();
2726            }
2727            delete [] pivotSave;
2728          }
2729#endif
2730        } else {
2731          bad=true;
2732        }
2733      }
2734      // clean up
2735      vector.setNumElements(0);
2736      for (int i=0;i<2*numberPivotsStored;i++) {
2737        array[--capacity]=0.0;
2738      }
2739      if (bad) {
2740        // switch off
2741#if ABC_NORMAL_DEBUG>0
2742        printf("switching off early factorization(in statusOf)\n");
2743#endif
2744#ifndef KEEP_TRYING_EARLY
2745        numberEarly_=0;
2746        delete abcEarlyFactorization_;
2747        abcEarlyFactorization_=NULL;
2748#else
2749        numberEarly_&=0xffff0000;
2750#endif
2751        factorizationStatus = internalFactorize(1);
2752      }
2753    }
2754#else
2755    factorizationStatus = internalFactorize(type ? 1 : 2);
2756#endif
2757    if (factorizationStatus) {
2758      if (type == 1&&lastFlaggedIteration_!=1) {
2759        // use for real disaster
2760        lastFlaggedIteration_ = 1;
2761        // no - restore previous basis
2762        unflagVariables = false;
2763        // restore weights
2764        weightsSaved = 4;
2765        changeMade_++; // say something changed
2766        // Keep any flagged variables
2767        for (int i = 0; i < numberTotal_; i++) {
2768          if (flagged(i))
2769            internalStatusSaved_[i] |= 64; //say flagged
2770        }
2771        restoreGoodStatus(type);
2772        // get correct bounds on all variables
2773        resetFakeBounds(1);
2774        numberBlackMarks=10;
2775        if (sequenceOut_>=0) {
2776          // need to reject something
2777          char x = isColumn(sequenceOut_) ? 'C' : 'R';
2778          handler_->message(CLP_SIMPLEX_FLAG, messages_)
2779            << x << sequenceWithin(sequenceOut_)
2780            << CoinMessageEol;
2781#if ABC_NORMAL_DEBUG>3
2782          printf("flag d\n");
2783#endif
2784          setFlagged(sequenceOut_);
2785          lastBadIteration_ = numberIterations_; // say be more cautious
2786          // so can't happen immediately again
2787          sequenceOut_=-1;
2788        }
2789        //abcProgress_.clearBadTimes();
2790       
2791        // Go to safe
2792        abcFactorization_->pivotTolerance(0.99);
2793      } else {
2794#if ABC_NORMAL_DEBUG>0
2795        printf("Bad initial basis\n");
2796#endif
2797      }
2798      if (internalFactorize(1)) {
2799        // try once more
2800      if (internalFactorize(0)) 
2801        abort();
2802      }
2803#if PARTITION_ROW_COPY==1
2804      int iVector=getAvailableArray();
2805      abcMatrix_->sortUseful(usefulArray_[iVector]);
2806      setAvailableArray(iVector);
2807#endif
2808    } else {
2809      lastFlaggedIteration_ = 0;
2810    }
2811  }
2812  // get primal and dual solutions
2813  numberBlackMarks+=gutsOfSolution(3)/2;
2814  assert (type!=0||numberIterations_==baseIteration_); // may be wrong
2815  if (!type) {
2816    initialSumInfeasibilities_=sumPrimalInfeasibilities_;
2817    initialNumberInfeasibilities_=numberPrimalInfeasibilities_;
2818    CoinAbcMemcpy(solutionSaved_+numberTotal_,abcSolution_,numberTotal_);
2819    CoinAbcMemcpy(internalStatusSaved_+numberTotal_,internalStatus_,numberTotal_);
2820  }
2821  double savePrimalInfeasibilities=sumPrimalInfeasibilities_;
2822  if ((moreSpecialOptions_&16384)!=0&&(moreSpecialOptions_&8192)==0) {
2823    if (numberIterations_==baseIteration_) {
2824      // If primal feasible and looks suited to primal go to primal
2825      if (!sumPrimalInfeasibilities_&&sumDualInfeasibilities_>1.0&&
2826          numberRows_*4>numberColumns_) {
2827        // do primal
2828        // mark so can perturb
2829        initialSumInfeasibilities_=-1.0;
2830        problemStatus_=10;
2831        return;
2832      }
2833    } else if ((numberIterations_>1000&&
2834               numberIterations_*10>numberRows_&&numberIterations_*4<numberRows_)||(sumPrimalInfeasibilities_+sumDualInfeasibilities_>
2835                                                                                    1.0e10*(initialSumInfeasibilities_+1000.0)&&numberIterations_>10000&&sumDualInfeasibilities_>1.0)) 
2836 {
2837      if (sumPrimalInfeasibilities_+1.0e2*sumDualInfeasibilities_>
2838          2000.0*(initialSumInfeasibilities_+100.0)) {
2839#if ABC_NORMAL_DEBUG>0
2840        printf ("sump %g sumd %g initial %g\n",sumPrimalInfeasibilities_,sumDualInfeasibilities_,
2841                initialSumInfeasibilities_);
2842#endif
2843        // but only if average gone up as well
2844        if (sumPrimalInfeasibilities_*initialNumberInfeasibilities_>
2845            100.0*(initialSumInfeasibilities_+10.0)*numberPrimalInfeasibilities_) {
2846#if ABC_NORMAL_DEBUG>0
2847          printf("returning to do primal\n");
2848#endif
2849          // do primal
2850          problemStatus_=10;
2851          CoinAbcMemcpy(abcSolution_,solutionSaved_+numberTotal_,numberTotal_);
2852          CoinAbcMemcpy(internalStatus_,internalStatusSaved_+numberTotal_,numberTotal_);
2853          // mark so can perturb
2854          initialSumInfeasibilities_=-1.0;
2855          return;
2856        }
2857      }
2858    }
2859  } else if ((sumPrimalInfeasibilities_+sumDualInfeasibilities_>
2860             1.0e6*(initialSumInfeasibilities_+1000.0)&&(moreSpecialOptions_&8192)==0&&
2861             numberIterations_>1000&&
2862              numberIterations_*10>numberRows_&&numberIterations_*4<numberRows_)||
2863             (sumPrimalInfeasibilities_+sumDualInfeasibilities_>
2864              1.0e10*(initialSumInfeasibilities_+1000.0)&&(moreSpecialOptions_&8192)==0&&
2865              numberIterations_>10000&&sumDualInfeasibilities_>1.0)) {
2866    // do primal
2867    problemStatus_=10;
2868    CoinAbcMemcpy(abcSolution_,solutionSaved_+numberTotal_,numberTotal_);
2869    CoinAbcMemcpy(internalStatus_,internalStatusSaved_+numberTotal_,numberTotal_);
2870    // mark so can perturb
2871    initialSumInfeasibilities_=-1.0;
2872    return;
2873  }
2874  //double saveDualInfeasibilities=sumDualInfeasibilities_;
2875  // 0 don't ignore, 1 ignore errors, 2 values pass
2876  int ignoreStuff=(firstFree_<0&&lastFirstFree_<0) ? 0 : 1;
2877  if ((stateOfProblem_&VALUES_PASS)!=0)
2878    ignoreStuff=2;
2879  {
2880    double thisObj = rawObjectiveValue();
2881    double lastObj = abcProgress_.lastObjective(0);
2882    if ((lastObj > thisObj + 1.0e-3 * CoinMax(fabs(thisObj), fabs(lastObj)) + 1.0e-4||
2883         (sumPrimalInfeasibilities_>10.0*lastSumPrimalInfeasibilities+1.0e3&&
2884          sumDualInfeasibilities_>10.0*CoinMax(lastSumDualInfeasibilities,1.0)))
2885        &&!ignoreStuff) {
2886#if ABC_NORMAL_DEBUG>0
2887      printf("bad - objective backwards %g %g errors %g %g suminf %g %g\n",lastObj,thisObj,largestDualError_,largestPrimalError_,sumDualInfeasibilities_,sumPrimalInfeasibilities_);
2888#endif
2889      numberBlackMarks++;
2890      if (largestDualError_>1.0e-3||largestPrimalError_>1.0e-3||
2891          (sumDualInfeasibilities_>10.0*CoinMax(lastSumDualInfeasibilities,1.0)
2892           &&firstFree_<0)) {
2893#if ABC_NORMAL_DEBUG>0
2894        printf("backtracking %d iterations - dual error %g primal error %g - sumDualInf %g\n",
2895               numberPivots,
2896               largestDualError_,largestPrimalError_,sumDualInfeasibilities_);
2897#endif
2898        // restore previous basis
2899        unflagVariables = false;
2900        // restore weights
2901        weightsSaved = 4;
2902        changeMade_++; // say something changed
2903        // Keep any flagged variables
2904        for (int i = 0; i < numberTotal_; i++) {
2905          if (flagged(i))
2906            internalStatusSaved_[i] |= 64; //say flagged
2907        }
2908        restoreGoodStatus(type);
2909        // get correct bounds on all variables
2910        resetFakeBounds(1);
2911        numberBlackMarks+=10;
2912        if (sequenceOut_>=0) {
2913          // need to reject something
2914#if ABC_NORMAL_DEBUG>0
2915          if (flagged(sequenceOut_))
2916            printf("BAD %x already flagged\n",sequenceOut_);
2917#endif
2918          char x = isColumn(sequenceOut_) ? 'C' : 'R';
2919          handler_->message(CLP_SIMPLEX_FLAG, messages_)
2920            << x << sequenceWithin(sequenceOut_)
2921            << CoinMessageEol;
2922#if ABC_NORMAL_DEBUG>3
2923          printf("flag d\n");
2924#endif
2925          setFlagged(sequenceOut_);
2926          lastBadIteration_ = numberIterations_; // say be more cautious
2927#if ABC_NORMAL_DEBUG>0
2928        } else {
2929          printf("backtracking - no valid sequence out\n");
2930#endif
2931        }
2932        // so can't happen immediately again
2933        sequenceOut_=-1;
2934        //abcProgress_.clearBadTimes();
2935        // Go to safe
2936        //abcFactorization_->pivotTolerance(0.99);
2937        if (internalFactorize(1)) {
2938          abort();
2939        }
2940#if PARTITION_ROW_COPY==1
2941        int iVector=getAvailableArray();
2942        abcMatrix_->sortUseful(usefulArray_[iVector]);
2943        setAvailableArray(iVector);
2944#endif
2945        // get primal and dual solutions
2946        numberBlackMarks+=gutsOfSolution(3);
2947      }
2948    }
2949  }
2950  //#define PAN
2951#ifdef PAN
2952    if (!type&&(specialOptions_&COIN_CBC_USING_CLP)==0&&firstFree_<0) {
2953    // first time - set large bad ones superbasic
2954    int whichArray=getAvailableArray();
2955    double * COIN_RESTRICT work=usefulArray_[whichArray].denseVector();
2956    int number=0;
2957    int * COIN_RESTRICT which=usefulArray_[whichArray].getIndices();
2958    const double * COIN_RESTRICT reducedCost=abcDj_;
2959    const double * COIN_RESTRICT lower=abcLower_;
2960    const double * COIN_RESTRICT upper=abcUpper_;
2961#if PRINT_PAN
2962    int nAtUb=0;
2963#endif
2964    for (int i=0;i<numberTotal_;i++) {
2965      Status status = getInternalStatus(i);
2966      double djValue=reducedCost[i];
2967      double badValue=0.0;
2968      if (status==atLowerBound&&djValue<-currentDualTolerance_) {
2969        double gap=upper[i]-lower[i];
2970        if (gap>100000000.0)
2971          badValue = -djValue*CoinMin(gap,1.0e10);
2972      } else if (status==atUpperBound&&djValue>currentDualTolerance_) {
2973        double gap=upper[i]-lower[i];
2974        if (gap>100000000.0) {
2975          badValue = djValue*CoinMin(gap,1.0e10);
2976#if PRINT_PAN
2977          nAtUb++;
2978#endif
2979        }
2980      }
2981      if (badValue) {
2982        work[number]=badValue;
2983        which[number++]=i;
2984      }
2985    }
2986    // for now don't do if primal feasible (probably should go to primal)
2987    if (!numberPrimalInfeasibilities_) {
2988      number=0;
2989#if ABC_NORMAL_DEBUG>1
2990      printf("XXXX doesn't work if primal feasible - also fix switch to Clp primal\n");
2991#endif
2992    }
2993    if (number&&ignoreStuff!=2) {
2994      CoinSort_2(work,work+number,which);
2995      int numberToDo=CoinMin(number,numberRows_/10); // ??
2996#if PRINT_PAN>1
2997      CoinSort_2(which,which+numberToDo,work);
2998      int n=0;
2999#endif
3000#if PRINT_PAN
3001      printf("Panning %d variables (out of %d at lb and %d at ub)\n",numberToDo,number-nAtUb,nAtUb);
3002#endif
3003      firstFree_=numberTotal_;
3004      for (int i=0;i<numberToDo;i++) {
3005        int iSequence=which[i];
3006#if PRINT_PAN>1
3007        if (!n)
3008          printf("Pan ");
3009        printf("%d ",iSequence);
3010        n++;
3011        if (n==10) {
3012          n=0;
3013          printf("\n");
3014        }
3015#endif
3016        setInternalStatus(iSequence,superBasic);
3017        firstFree_=CoinMin(firstFree_,iSequence);
3018      }
3019#if PRINT_PAN>1
3020      if (n)
3021        printf("\n");
3022#endif
3023      ordinaryVariables_=0;
3024      CoinAbcMemset0(work,number);
3025      stateOfProblem_ |= FAKE_SUPERBASIC;
3026    }
3027    setAvailableArray(whichArray);
3028  }
3029#endif
3030  if (!sumOfRelaxedPrimalInfeasibilities_&&sumDualInfeasibilities_&&type&&saveType<9
3031      &&ignoreStuff!=2) {
3032    // say been feasible
3033    abcState_ |= CLP_ABC_BEEN_FEASIBLE;
3034    int numberFake=numberAtFakeBound();
3035    if (!numberFake) {
3036      makeNonFreeVariablesDualFeasible();
3037      numberDualInfeasibilities_=0;
3038      sumDualInfeasibilities_=0.0;
3039      // just primal
3040      gutsOfSolution(2);
3041    }
3042  }
3043  bool computeNewSumInfeasibilities=numberIterations_==baseIteration_;
3044  if ((!type||(firstFree_<0&&lastFirstFree_>=0))&&ignoreStuff!=2) {
3045#ifdef PAN
3046    if ((stateOfProblem_&FAKE_SUPERBASIC)!=0&&firstFree_<0&&lastFirstFree_>=0&&type) {
3047      // clean up
3048#if PRINT_PAN
3049      printf("Pan switch off after %d iterations\n",numberIterations_);
3050#endif
3051      computeNewSumInfeasibilities=true;
3052      ordinaryVariables_=1;
3053      const double * COIN_RESTRICT reducedCost=abcDj_;
3054      const double * COIN_RESTRICT lower=abcLower_;
3055      const double * COIN_RESTRICT upper=abcUpper_;
3056      for (int i=0;i<numberTotal_;i++) {
3057        Status status = getInternalStatus(i);
3058        if (status==superBasic) {
3059          double djValue=reducedCost[i];
3060          double value=abcSolution_[i];
3061          if (value<lower[i]+primalTolerance_) {
3062            setInternalStatus(i,atLowerBound);
3063#if ABC_NORMAL_DEBUG>0
3064            if(djValue<-dualTolerance_) 
3065              printf("Odd at lb %d dj %g\n",i,djValue);
3066#endif
3067          } else if (value>upper[i]-primalTolerance_) {
3068            setInternalStatus(i,atUpperBound);
3069#if ABC_NORMAL_DEBUG>0
3070            if(djValue>dualTolerance_) 
3071              printf("Odd at ub %d dj %g\n",i,djValue);
3072#endif
3073          } else {
3074            assert (status!=superBasic);
3075          }
3076        }
3077        if (status==isFree)
3078          ordinaryVariables_=0;
3079      }
3080      stateOfProblem_ &= ~FAKE_SUPERBASIC;
3081    }
3082#endif
3083    // set up perturbation and dualBound_
3084    // make dual feasible
3085    if (firstFree_<0) {
3086      if (numberFlagged_) {
3087        numberFlagged_=0;
3088        for (int iRow = 0; iRow < numberRows_; iRow++) {
3089          int iPivot = abcPivotVariable_[iRow];
3090          if (flagged(iPivot)) {
3091            clearFlagged(iPivot);
3092          }
3093        }
3094      }
3095      int numberChanged=resetFakeBounds(0);
3096#if ABC_NORMAL_DEBUG>0
3097      if (numberChanged>0)
3098        printf("Re-setting %d fake bounds %d gone dual infeasible - can expect backward objective\n",
3099               numberChanged,numberDualInfeasibilities_);
3100#endif
3101      if (numberDualInfeasibilities_||numberChanged!=-1) {
3102        makeNonFreeVariablesDualFeasible();
3103        numberDualInfeasibilities_=0;
3104        sumDualInfeasibilities_=0.0;
3105        sumOfRelaxedDualInfeasibilities_=0.0;
3106        // just primal
3107        gutsOfSolution(2);
3108        if (!numberIterations_ && sumPrimalInfeasibilities_ >
3109            1.0e5*(savePrimalInfeasibilities+1.0e3) &&
3110            (moreSpecialOptions_ & (256|8192)) == 0) {
3111          // Use primal
3112          problemStatus_=10;
3113          // mark so can perturb
3114          initialSumInfeasibilities_=-1.0;
3115          CoinAbcMemcpy(abcSolution_,solutionSaved_,numberTotal_);
3116          CoinAbcMemcpy(internalStatus_,internalStatusSaved_,numberTotal_);
3117          return;
3118        }
3119        // redo
3120        if(computeNewSumInfeasibilities)
3121          initialSumInfeasibilities_=sumPrimalInfeasibilities_;
3122        //computeObjective();
3123      } 
3124    }
3125    //bounceTolerances(type);
3126  } else if (numberIterations_>numberRows_&&stateDualColumn_==-2) {
3127    bounceTolerances(4);
3128  }
3129  // If bad accuracy treat as singular
3130  if ((largestPrimalError_ > 1.0e15 || largestDualError_ > 1.0e15) && numberIterations_>1000) {
3131    problemStatus_=10; //abort();
3132#if ABC_NORMAL_DEBUG>0
3133    printf("Big problems - errors %g %g try primal\n",largestPrimalError_,largestDualError_);
3134#endif
3135  } else if (goodAccuracy()&&numberIterations_>lastBadIteration_+200) {
3136    // Can reduce tolerance
3137    double newTolerance = CoinMax(0.995 * abcFactorization_->pivotTolerance(), saveData_.pivotTolerance_);
3138    if (!type)
3139      newTolerance = saveData_.pivotTolerance_;
3140    if (newTolerance>abcFactorization_->minimumPivotTolerance())
3141      abcFactorization_->pivotTolerance(newTolerance);
3142  }
3143  bestObjectiveValue_ = CoinMax(bestObjectiveValue_,
3144                                rawObjectiveValue_ - bestPossibleImprovement_);
3145  // Double check infeasibility if no action
3146  if (abcProgress_.lastIterationNumber(0) == numberIterations_) {
3147    if (abcDualRowPivot_->looksOptimal()) {
3148      // say been feasible
3149      abcState_ |= CLP_ABC_BEEN_FEASIBLE;
3150      numberPrimalInfeasibilities_ = 0;
3151      sumPrimalInfeasibilities_ = 0.0;
3152    }
3153  } else {
3154    double thisObj = rawObjectiveValue_ - bestPossibleImprovement_;
3155#ifdef CLP_INVESTIGATE
3156    assert (bestPossibleImprovement_ > -1000.0 && rawObjectiveValue_ > -1.0e100);
3157    if (bestPossibleImprovement_)
3158      printf("obj %g add in %g -> %g\n", rawObjectiveValue_, bestPossibleImprovement_,
3159             thisObj);
3160#endif
3161    double lastObj = abcProgress_.lastObjective(0);
3162#ifndef NDEBUG
3163#if ABC_NORMAL_DEBUG>3
3164    resetFakeBounds(-1);
3165#endif
3166#endif
3167    if (lastObj < thisObj - 1.0e-5 * CoinMax(fabs(thisObj), fabs(lastObj)) - 1.0e-3 && saveType<9) {
3168      numberTimesOptimal_ = 0;
3169    }
3170  }
3171  // Up tolerance if looks a bit odd
3172  if (numberIterations_ > CoinMax(1000, numberRows_ >> 4) && (specialOptions_ & 64) != 0) {
3173    if (sumPrimalInfeasibilities_ && sumPrimalInfeasibilities_ < 1.0e5) {
3174      int backIteration = abcProgress_.lastIterationNumber(CLP_PROGRESS - 1);
3175      if (backIteration > 0 && numberIterations_ - backIteration < 9 * CLP_PROGRESS) {
3176        if (abcFactorization_->pivotTolerance() < 0.9) {
3177          // up tolerance
3178          abcFactorization_->pivotTolerance(CoinMin(abcFactorization_->pivotTolerance() * 1.05 + 0.02, 0.91));
3179          //printf("tol now %g\n",abcFactorization_->pivotTolerance());
3180          abcProgress_.clearIterationNumbers();
3181        }
3182      }
3183    }
3184  }
3185  // Check if looping
3186  int loop;
3187  if (type != 2)
3188    loop = abcProgress_.looping();
3189  else
3190    loop = -1;
3191  if (abcProgress_.reallyBadTimes() > 10) {
3192    problemStatus_ = 10; // instead - try other algorithm
3193#if ABC_NORMAL_DEBUG>3
3194    printf("returning at %d\n", __LINE__);
3195#endif
3196  }
3197  if (loop >= 0) {
3198    problemStatus_ = loop; //exit if in loop
3199    if (!problemStatus_) {
3200      // declaring victory
3201      // say been feasible
3202      abcState_ |= CLP_ABC_BEEN_FEASIBLE;
3203      numberPrimalInfeasibilities_ = 0;
3204      sumPrimalInfeasibilities_ = 0.0;
3205    } else {
3206      problemStatus_ = 10; // instead - try other algorithm
3207#if ABC_NORMAL_DEBUG>3
3208      printf("returning at %d\n", __LINE__);
3209#endif
3210    }
3211    return;
3212  }
3213  // really for free variables in
3214  if((progressFlag_ & 2) != 0) {
3215    //situationChanged = 2;
3216  }
3217  progressFlag_ &= (~3); //reset progress flag
3218  // mark as having gone optimal if looks like it
3219  if (!numberPrimalInfeasibilities_&&
3220      !numberDualInfeasibilities_)
3221    progressFlag_ |= 8;
3222  if (handler_->detail(CLP_SIMPLEX_STATUS, messages_) < 100) {
3223    handler_->message(CLP_SIMPLEX_STATUS, messages_)
3224      << numberIterations_ << objectiveValue();
3225    handler_->printing(sumPrimalInfeasibilities_ > 0.0)
3226      << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
3227    handler_->printing(sumDualInfeasibilities_ > 0.0)
3228      << sumDualInfeasibilities_ << numberDualInfeasibilities_;
3229    handler_->printing(numberDualInfeasibilitiesWithoutFree_
3230                       < numberDualInfeasibilities_)
3231      << numberDualInfeasibilitiesWithoutFree_;
3232    handler_->message() << CoinMessageEol;
3233  }
3234  double approximateObjective = rawObjectiveValue_;
3235  //realDualInfeasibilities = sumDualInfeasibilities_;
3236  // If we need to carry on cleaning variables
3237  if (!numberPrimalInfeasibilities_ && (specialOptions_ & 1024) != 0 && CLEAN_FIXED) {
3238    for (int iRow = 0; iRow < numberRows_; iRow++) {
3239      int iPivot = abcPivotVariable_[iRow];
3240      if (!flagged(iPivot) && pivoted(iPivot)) {
3241        // carry on
3242        numberPrimalInfeasibilities_ = -1;
3243        sumOfRelaxedPrimalInfeasibilities_ = 1.0;
3244        sumPrimalInfeasibilities_ = 1.0;
3245        break;
3246      }
3247    }
3248  }
3249  /* If we are primal feasible and any dual infeasibilities are on
3250     free variables then it is better to go to primal */
3251  if (!numberPrimalInfeasibilities_ && !numberDualInfeasibilitiesWithoutFree_ &&
3252      numberDualInfeasibilities_) {
3253    // say been feasible
3254    abcState_ |= CLP_ABC_BEEN_FEASIBLE;
3255    if (type&&ignoreStuff!=2)
3256      problemStatus_ = 10;
3257    else
3258      numberPrimalInfeasibilities_=1;
3259  }
3260  // check optimal
3261  // give code benefit of doubt
3262  if (sumOfRelaxedDualInfeasibilities_ == 0.0 &&
3263      sumOfRelaxedPrimalInfeasibilities_ == 0.0&&ignoreStuff!=2) {
3264    // say been feasible
3265    abcState_ |= CLP_ABC_BEEN_FEASIBLE;
3266    // say optimal (with these bounds etc)
3267    numberDualInfeasibilities_ = 0;
3268    sumDualInfeasibilities_ = 0.0;
3269    numberPrimalInfeasibilities_ = 0;
3270    sumPrimalInfeasibilities_ = 0.0;
3271  } else if (sumOfRelaxedDualInfeasibilities_>1.0e-3&&firstFree_<0&&ignoreStuff!=2
3272             &&numberIterations_>lastCleaned_+10) {
3273#if ABC_NORMAL_DEBUG>0
3274    printf("Making dual feasible\n");
3275#endif
3276    makeNonFreeVariablesDualFeasible(true);
3277    lastCleaned_=numberIterations_;
3278    numberDualInfeasibilities_=0;
3279    sumDualInfeasibilities_=0.0;
3280    // just primal
3281    gutsOfSolution(2);
3282  }
3283  // see if cutoff reached
3284  double limit = 0.0;
3285  getDblParam(ClpDualObjectiveLimit, limit);
3286  if (primalFeasible()&&ignoreStuff!=2) {
3287    // may be optimal - or may be bounds are wrong
3288    int numberChangedBounds = changeBounds(0, changeCost);
3289    if (numberChangedBounds <= 0 && !numberDualInfeasibilities_) {
3290      //looks optimal - do we need to reset tolerance
3291      if (lastCleaned_ < numberIterations_ && /*numberTimesOptimal_*/ stateDualColumn_ < 4 &&
3292          (numberChanged_ || (specialOptions_ & 4096) == 0)) {
3293#if CLP_CAN_HAVE_ZERO_OBJ
3294        if ((specialOptions_&2097152)==0) {
3295#endif
3296          numberTimesOptimal_++;
3297          if (stateDualColumn_==-2) {
3298#if ABC_NORMAL_DEBUG>0
3299            printf("Going straight from state -2 to 0\n");
3300#endif
3301            stateDualColumn_=-1;
3302          }
3303          changeMade_++; // say something changed
3304          if (numberTimesOptimal_ == 1) {
3305            if (fabs(currentDualTolerance_-dualTolerance_)>1.0e-12) {
3306#if ABC_NORMAL_DEBUG>0
3307              printf("changing dual tolerance from %g to %g (numberTimesOptimal_==%d)\n",
3308                     currentDualTolerance_,dualTolerance_,numberTimesOptimal_);
3309#endif
3310              currentDualTolerance_ = dualTolerance_;
3311            }
3312          } else {
3313            if (numberTimesOptimal_ == 2) {
3314              // better to have small tolerance even if slower
3315              abcFactorization_->zeroTolerance(CoinMin(abcFactorization_->zeroTolerance(), 1.0e-15));
3316            }
3317#if ABC_NORMAL_DEBUG>0
3318            printf("changing dual tolerance from %g to %g (numberTimesOptimal_==%d)\n",
3319                   currentDualTolerance_,dualTolerance_*pow(2.0,numberTimesOptimal_-1),numberTimesOptimal_);
3320#endif
3321            currentDualTolerance_ = dualTolerance_ * pow(2.0, numberTimesOptimal_ - 1);
3322          }
3323          if (numberChanged_>0) {
3324            handler_->message(CLP_DUAL_ORIGINAL, messages_)
3325              << CoinMessageEol;
3326            //realDualInfeasibilities = sumDualInfeasibilities_;
3327          } else if (numberChangedBounds<-1) {
3328            // bounds were flipped
3329            //gutsOfSolution(2);
3330          }
3331          if (stateDualColumn_>=0)
3332            stateDualColumn_++;
3333#ifndef TRY_ABC_GUS
3334          int returnCode=bounceTolerances(1);
3335          if (returnCode)
3336#endif // always go to primal
3337            problemStatus_=10;
3338          lastCleaned_ = numberIterations_;
3339#if CLP_CAN_HAVE_ZERO_OBJ
3340        } else {
3341          // no cost - skip checks
3342          problemStatus_=0;
3343        }
3344#endif
3345      } else {
3346        problemStatus_ = 0; // optimal
3347        if (lastCleaned_ < numberIterations_ && numberChanged_) {
3348          handler_->message(CLP_SIMPLEX_GIVINGUP, messages_)
3349            << CoinMessageEol;
3350        }
3351      }
3352    } else if (numberChangedBounds<-1) {
3353      // some flipped
3354      bounceTolerances(2);
3355#if ABC_NORMAL_DEBUG>0
3356      printf("numberChangedBounds %d numberAtFakeBound %d at line %d in file %s\n",
3357             numberChangedBounds,numberAtFakeBound(),__LINE__,__FILE__);
3358#endif
3359    } else if (numberAtFakeBound()) {
3360      handler_->message(CLP_DUAL_BOUNDS, messages_)
3361        << currentDualBound_
3362        << CoinMessageEol;
3363#if ABC_NORMAL_DEBUG>0
3364      printf("changing bounds\n");
3365      printf("numberChangedBounds %d numberAtFakeBound %d at line %d in file %s\n",
3366             numberChangedBounds,numberAtFakeBound(),__LINE__,__FILE__);
3367#endif
3368      if (!numberDualInfeasibilities_&&!numberPrimalInfeasibilities_) {
3369        // check unbounded
3370        // find a variable with bad dj
3371        int iChosen = -1;
3372        if ((specialOptions_ & 0x03000000) == 0) {
3373          double largest = 100.0 * primalTolerance_;
3374          for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
3375            double djValue = abcDj_[iSequence];
3376            double originalLo = lowerSaved_[iSequence];
3377            double originalUp = upperSaved_[iSequence];
3378            if (fabs(djValue) > fabs(largest)) {
3379              if (getInternalStatus(iSequence) != basic) {
3380                if (djValue > 0 && originalLo < -1.0e20) {
3381                  if (djValue > fabs(largest)) {
3382                    largest = djValue;
3383                    iChosen = iSequence;
3384                  }
3385                } else if (djValue < 0 && originalUp > 1.0e20) {
3386                  if (-djValue > fabs(largest)) {
3387                    largest = djValue;
3388                    iChosen = iSequence;
3389                  }
3390                }
3391              }
3392            }
3393          }
3394        }
3395        if (iChosen >= 0) {
3396          int iVector=getAvailableArray();
3397          unpack(usefulArray_[iVector],iChosen);
3398          //double changeCost;
3399          problemStatus_ = checkUnbounded(usefulArray_[iVector], 0.0/*changeCost*/);
3400          usefulArray_[iVector].clear();
3401          setAvailableArray(iVector);
3402        }
3403      }
3404      gutsOfSolution(3);
3405      //bounceTolerances(2);
3406      //realDualInfeasibilities = sumDualInfeasibilities_;
3407    }
3408    // unflag all variables (we may want to wait a bit?)
3409    if (unflagVariables) {
3410      int numberFlagged = numberFlagged_;
3411      numberFlagged_=0;
3412      for (int iRow = 0; iRow < numberRows_; iRow++) {
3413        int iPivot = abcPivotVariable_[iRow];
3414        if (flagged(iPivot)) {
3415          clearFlagged(iPivot);
3416        }
3417      }
3418#if ABC_NORMAL_DEBUG>3
3419      if (numberFlagged) {
3420        printf("unflagging %d variables - tentativeStatus %d probStat %d ninf %d nopt %d\n", numberFlagged, tentativeStatus,
3421               problemStatus_, numberPrimalInfeasibilities_,
3422               numberTimesOptimal_);
3423      }
3424#endif
3425      unflagVariables = numberFlagged > 0;
3426      if (numberFlagged && !numberPivots) {
3427        /* looks like trouble as we have not done any iterations.
3428           Try changing pivot tolerance then give it a few goes and give up */
3429        if (abcFactorization_->pivotTolerance() < 0.9) {
3430          abcFactorization_->pivotTolerance(0.99);
3431          problemStatus_ = -1;
3432        } else if (numberTimesOptimal_ < 3) {
3433          numberTimesOptimal_++;
3434          problemStatus_ = -1;
3435        } else {
3436          unflagVariables = false;
3437          //secondaryStatus_ = 1; // and say probably infeasible
3438          if ((moreSpecialOptions_ & 256) == 0||(abcState_&CLP_ABC_BEEN_FEASIBLE)!=0) {
3439            // try primal
3440            problemStatus_ = 10;
3441          } else {
3442            // almost certainly infeasible
3443            problemStatus_ = 1;
3444          }
3445#if ABC_NORMAL_DEBUG>3
3446          printf("returning at %d\n", __LINE__);
3447#endif
3448        }
3449      }
3450    }
3451  }
3452  if (problemStatus_ < 0) {
3453    saveGoodStatus();
3454    if (weightsSaved) {
3455      // restore weights (if saved) - also recompute infeasibility list
3456      abcDualRowPivot_->saveWeights(this, weightsSaved);
3457    } else {
3458      abcDualRowPivot_->recomputeInfeasibilities();
3459    }
3460  }
3461  // see if cutoff reached
3462  checkCutoff(false);
3463  // If we are in trouble and in branch and bound give up
3464  if ((specialOptions_ & 1024) != 0) {
3465    int looksBad = 0;
3466    if (largestPrimalError_ * largestDualError_ > 1.0e2) {
3467      looksBad = 1;
3468    } else if (largestPrimalError_ > 1.0e-2
3469               && rawObjectiveValue_ > CoinMin(1.0e15, 1.0e3 * limit)) {
3470      looksBad = 2;
3471    }
3472    if (looksBad) {
3473      if (abcFactorization_->pivotTolerance() < 0.9) {
3474        // up tolerance
3475        abcFactorization_->pivotTolerance(CoinMin(abcFactorization_->pivotTolerance() * 1.05 + 0.02, 0.91));
3476      } else if (numberIterations_ > 10000) {
3477#if ABC_NORMAL_DEBUG>0
3478        if (handler_->logLevel() > 2)
3479          printf("bad dual - saying infeasible %d\n", looksBad);
3480#endif
3481        problemStatus_ = 1;
3482        secondaryStatus_ = 1; // and say was on cutoff
3483      } else if (largestPrimalError_ > 1.0e5) {
3484#if ABC_NORMAL_DEBUG>0
3485        if (handler_->logLevel() > 2)
3486          printf("bad dual - going to primal %d %g\n", looksBad, largestPrimalError_);
3487#endif
3488        allSlackBasis();
3489        problemStatus_ = 10;
3490      }
3491    }
3492  }
3493  //if (problemStatus_ < 0 && !changeMade_) {
3494  //problemStatus_ = 4; // unknown
3495  //}
3496  lastGoodIteration_ = numberIterations_;
3497  if (numberIterations_ > lastBadIteration_ + 200) {
3498    moreSpecialOptions_ &= ~16; // clear check accuracy flag
3499    if (numberFlagged_) {
3500      numberFlagged_=0;
3501      for (int iRow = 0; iRow < numberRows_; iRow++) {
3502        int iPivot = abcPivotVariable_[iRow];
3503        if (flagged(iPivot)) {
3504          clearFlagged(iPivot);
3505        }
3506      }
3507    }
3508  }
3509  if (problemStatus_ < 0&&ignoreStuff!=2) {
3510    //sumDualInfeasibilities_ = realDualInfeasibilities; // back to say be careful
3511    if (sumDualInfeasibilities_) {
3512      numberDualInfeasibilities_ = 1;
3513    } else if (!numberPrimalInfeasibilities_) {
3514      // Looks good
3515      problemStatus_=0;
3516    }
3517  }
3518  double thisObj = abcProgress_.lastObjective(0);
3519  double lastObj = abcProgress_.lastObjective(1);
3520  if (lastObj > thisObj + 1.0e-4 * CoinMax(fabs(thisObj), fabs(lastObj)) + 1.0e-4&&
3521      lastFirstFree_<0) {
3522    //printf("BAD - objective backwards\n");
3523    //bounceTolerances(3);
3524    numberBlackMarks+=3;
3525  }
3526  if (numberBlackMarks>0&&numberIterations_>baseIteration_) {
3527    // be very careful
3528    int maxFactor = abcFactorization_->maximumPivots();
3529    if (maxFactor > 10) {
3530      if ((stateOfProblem_&PESSIMISTIC)==0||true) {
3531        if (largestDualError_>10.0*CoinMax(lastDualError_,1.0e-6)
3532            ||largestPrimalError_>10.0*CoinMax(lastPrimalError_,1.0e-6))
3533          maxFactor = CoinMin(maxFactor,20);
3534        forceFactorization_ = CoinMax(1, (maxFactor >> 1));
3535      } else if (numberBlackMarks>2) {
3536        forceFactorization_ = CoinMax(1, (forceFactorization_ >> 1));
3537      }
3538    }
3539    stateOfProblem_ |= PESSIMISTIC;
3540    if (numberBlackMarks>=10)
3541      forceFactorization_ = 1;
3542  }
3543  lastPrimalError_=largestPrimalError_;
3544  lastDualError_=largestDualError_;
3545  lastFirstFree_=firstFree_;
3546  if (ignoreStuff==2) {
3547    // in values pass
3548    lastFirstFree_=-1;
3549    // firstFree_=-1;
3550  }
3551  if (alphaAccuracy_ > 0.0)
3552    alphaAccuracy_ = 1.0;
3553  // If we are stopping - use plausible objective
3554  // Maybe only in fast dual
3555  if (problemStatus_ > 2)
3556    rawObjectiveValue_ = approximateObjective;
3557  if (problemStatus_ == 1 && (progressFlag_&8) != 0 &&
3558      fabs(rawObjectiveValue_) > 1.0e10 )
3559    problemStatus_ = 10; // infeasible - but has looked feasible
3560  checkMoveBack(true);
3561#if 0
3562  {
3563    double * array = usefulArray_[arrayForTableauRow_].denseVector();
3564    for (int iPivot = 0; iPivot < numberRows_; iPivot++) {
3565      int iSequence = abcPivotVariable_[iPivot];
3566      unpack(usefulArray_[arrayForTableauRow_],iSequence);
3567      abcFactorization_->updateColumn(usefulArray_[arrayForTableauRow_]);
3568      assert (fabs(array[iPivot] - 1.0) < 1.0e-4);
3569      array[iPivot] = 0.0;
3570      for (int i = 0; i < numberRows_; i++)
3571        assert (fabs(array[i]) < 1.0e-4);
3572      usefulArray_[arrayForTableauRow_].clear();
3573    }
3574  }
3575#endif
3576#if 0
3577  {
3578    // save status, cost and solution
3579    unsigned char * saveStatus = CoinCopyOfArray(internalStatus_,numberTotal_);
3580    double * saveSolution = CoinCopyOfArray(abcSolution_,numberTotal_);
3581    double * saveCost = CoinCopyOfArray(abcCost_,numberTotal_);
3582    printf("TestA sumDual %g sumPrimal %g\n",sumDualInfeasibilities_,sumPrimalInfeasibilities_);
3583    if (perturbation_>100)
3584      CoinAbcMemcpy(abcCost_,costSaved_,numberTotal_);
3585    else
3586      CoinAbcMemcpy(abcCost_,abcPerturbation_,numberTotal_);
3587    moveToBasic();
3588    gutsOfSolution(3);
3589    for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
3590      bool changed=fabs(saveCost[iSequence]-abcCost_[iSequence])>1.0e-10;
3591      bool bad=false;
3592      int iStatus=internalStatus_[iSequence]&7;
3593      if (getInternalStatus(iSequence) != basic && !flagged(iSequence)) {
3594        // not basic
3595        double distanceUp = abcUpper_[iSequence] - abcSolution_[iSequence];
3596        double distanceDown = abcSolution_[iSequence] - abcLower_[iSequence];
3597        double value = abcDj_[iSequence];
3598        if (distanceUp > primalTolerance_) {
3599          // Check if "free"
3600          if (distanceDown <= primalTolerance_) {
3601            // should not be negative
3602            if (value < 0.0) {
3603              value = - value;
3604              if (value > currentDualTolerance_) {
3605                bad=true;
3606              }
3607            }
3608          } else {
3609            // free
3610            value=fabs(value);
3611            if (value > currentDualTolerance_) {
3612              bad=true;
3613            }
3614          }
3615        } else if (distanceDown > primalTolerance_) {
3616          // should not be positive
3617          if (value > 0.0) {
3618            if (value > currentDualTolerance_) {
3619              bad=true;
3620            }
3621          }
3622        }
3623      }
3624      if (changed||bad)
3625        printf("seq %d status %d current cost %g original %g dj %g %s\n",
3626               iSequence,iStatus,saveCost[iSequence],abcCost_[iSequence],
3627               abcDj_[iSequence],bad ? "(bad)" : "");
3628    }
3629    printf("TestB sumDual %g sumPrimal %g\n",sumDualInfeasibilities_,sumPrimalInfeasibilities_);
3630    if (numberDualInfeasibilities_) {
3631      makeNonFreeVariablesDualFeasible();
3632      moveToBasic();
3633      gutsOfSolution(3);
3634      printf("TestC sumDual %g sumPrimal %g\n",sumDualInfeasibilities_,sumPrimalInfeasibilities_);
3635    }
3636    CoinAbcMemcpy(internalStatus_,saveStatus,numberTotal_);
3637    CoinAbcMemcpy(abcSolution_,saveSolution,numberTotal_);
3638    CoinAbcMemcpy(abcCost_,saveCost,numberTotal_);
3639    delete [] saveStatus;
3640    delete [] saveSolution;
3641    delete [] saveCost;
3642    moveToBasic();
3643    gutsOfSolution(3);
3644    printf("TestD sumDual %g sumPrimal %g\n",sumDualInfeasibilities_,sumPrimalInfeasibilities_);
3645  }
3646#endif
3647#if 0
3648  {
3649    const double multiplier[] = { 1.0, -1.0}; // sign?????????????
3650    double dualT = - currentDualTolerance_;
3651    double largestCostDifference=0.1*currentDualTolerance_;
3652    int numberCostDifference=0;
3653    double largestBadDj=0.0;
3654    int numberBadDj=0;
3655    double sumBadDj=0.0;
3656    double sumCanMoveBad=0.0;
3657    int numberCanMoveBad=0;
3658    double sumCanMoveOthers=0.0;
3659    int numberCanMoveOthers=0;
3660    for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
3661      double costDifference=abcPerturbation_[iSequence]-abcCost_[iSequence];
3662      if (fabs(costDifference)>0.1*currentDualTolerance_) {
3663        if (fabs(costDifference)>largestCostDifference)
3664          largestCostDifference=fabs(costDifference);
3665        numberCostDifference++;
3666      }
3667      unsigned char iStatus = internalStatus_[iSequence] & 7;
3668      if (iStatus<4) {
3669        double mult = multiplier[iStatus];
3670        double dj = abcDj_[iSequence] * mult;
3671        if (dj < dualT) {
3672          sumBadDj+=dualT-dj;
3673          numberBadDj++;
3674          if (dualT-dj>largestBadDj)
3675            largestBadDj=dualT-dj;
3676          if (costDifference*mult>0.0) {
3677            sumCanMoveBad+=costDifference*mult;
3678            numberCanMoveBad++;
3679            //dj+=costDifference*mult;
3680            abcDj_[iSequence]+=costDifference;
3681            abcCost_[iSequence]=abcPerturbation_[iSequence];
3682          }
3683        } else if (costDifference*mult>0.0) {
3684          sumCanMoveOthers+=costDifference*mult;
3685          numberCanMoveOthers++;
3686          //dj+=costDifference*mult;
3687          abcDj_[iSequence]+=costDifference;
3688          abcCost_[iSequence]=abcPerturbation_[iSequence];
3689        }
3690      }
3691    }
3692    if (numberCostDifference+numberBadDj+numberCanMoveBad+numberCanMoveOthers) {
3693      printf("largest diff %g (%d), largest bad dj %g (%d - sum %g), can move bad %g (%d), good %g (%d) - %d at fake bound\n",largestCostDifference,numberCostDifference,largestBadDj,
3694             numberBadDj, sumBadDj,sumCanMoveBad=0.0,numberCanMoveBad,
3695             sumCanMoveOthers,numberCanMoveOthers,numberAtFakeBound());
3696    }
3697  }
3698#endif
3699}
3700// see if cutoff reached
3701bool
3702AbcSimplexDual::checkCutoff(bool computeObjective)
3703{
3704  double objValue=COIN_DBL_MAX;
3705  if (computeObjective)
3706    objValue = computeInternalObjectiveValue();
3707  else
3708    objValue = minimizationObjectiveValue();
3709  double limit = 0.0;
3710  getDblParam(ClpDualObjectiveLimit, limit);
3711  // Looks as if numberAtFakeBound() being computed too often
3712  if(fabs(limit) < 1.0e30 && objValue >
3713     limit && !numberAtFakeBound()) {
3714    bool looksInfeasible = !numberDualInfeasibilities_;
3715    if (objValue > limit + fabs(0.1 * limit) + 1.0e2 * sumDualInfeasibilities_ + 1.0e4 &&
3716        sumDualInfeasibilities_ < largestDualError_ && numberIterations_ > 0.5 * numberRows_ + 1000)
3717      looksInfeasible = true;
3718    if (looksInfeasible&&!computeObjective)
3719      objValue = computeInternalObjectiveValue();
3720    if (looksInfeasible) {
3721      // Even if not perturbed internal costs may have changed
3722      // be careful
3723      if(objValue > limit) {
3724        problemStatus_ = 1;
3725        secondaryStatus_ = 1; // and say was on cutoff
3726        return true;
3727      }
3728    }
3729  }
3730  return false;
3731}
3732#define CLEANUP_DJS 0
3733/* While updateDualsInDual sees what effect is of flip
3734   this does actual flipping.
3735   If change >0.0 then value in array >0.0 => from lower to upper
3736   returns 3 if skip this iteration and re-factorize
3737*/
3738int
3739AbcSimplexDual::flipBounds()
3740{
3741  CoinIndexedVector & array = usefulArray_[arrayForFlipBounds_];
3742  CoinIndexedVector & output = usefulArray_[arrayForFlipRhs_];
3743  output.clear();
3744  int numberFlipped=array.getNumElements();
3745  if (numberFlipped) {
3746#if CLP_CAN_HAVE_ZERO_OBJ>1
3747    if ((specialOptions_&2097152)==0) {
3748#endif
3749      // do actual flips
3750      // do in parallel once got spare space in factorization
3751#if ABC_PARALLEL==2
3752#endif
3753      int number=array.getNumElements();
3754      const int *  COIN_RESTRICT which=array.getIndices();
3755      double *  COIN_RESTRICT work = array.denseVector();
3756     
3757      double *  COIN_RESTRICT solution = abcSolution_;
3758      const double *  COIN_RESTRICT lower = abcLower_;
3759      const double *  COIN_RESTRICT upper = abcUpper_;
3760#if CLEANUP_DJS
3761      double *  COIN_RESTRICT cost = abcCost_;
3762      double *  COIN_RESTRICT dj = abcDj_;
3763      // see if we can clean up djs
3764      // may have perturbation
3765      const double * COIN_RESTRICT originalCost = abcPerturbation_;
3766#else
3767      const double *  COIN_RESTRICT cost = abcCost_;
3768      const double *  COIN_RESTRICT dj = abcDj_;
3769#endif
3770      const double multiplier[] = { 1.0, -1.0};
3771      for (int i = 0; i < number; i++) {
3772        double value=work[i]*theta_;
3773        work[i]=0.0;
3774        int iSequence = which[i];
3775        unsigned char iStatus = internalStatus_[iSequence]&3;
3776        assert ((internalStatus_[iSequence]&7)==0||
3777                (internalStatus_[iSequence]&7)==1);
3778        double mult = multiplier[iStatus];
3779        double djValue = dj[iSequence] * mult-value;
3780        //double perturbedTolerance = abcPerturbation_[iSequence];
3781        if (djValue<-currentDualTolerance_) {
3782          // might just happen - if so just skip
3783          assert (iSequence!=sequenceIn_);
3784          double movement=(upper[iSequence]-lower[iSequence])*mult;
3785          if (iStatus==0) {
3786            // at lower bound
3787            // to upper bound
3788            setInternalStatus(iSequence, atUpperBound);
3789            solution[iSequence] = upper[iSequence];
3790#if CLEANUP_DJS
3791            if (cost[iSequence]<originalCost[iSequence]) {
3792              double difference = CoinMax(cost[iSequence]-originalCost[iSequence],djValue);
3793              dj[iSequence]-=difference;
3794              cost[iSequence]-=difference;
3795            }
3796#endif
3797          } else {
3798            // at upper bound
3799            // to lower bound
3800            setInternalStatus(iSequence, atLowerBound);
3801            solution[iSequence] = lower[iSequence];
3802#if CLEANUP_DJS
3803            if (cost[iSequence]>originalCost[iSequence]) {
3804              double difference = CoinMin(cost[iSequence]-originalCost[iSequence],-djValue);
3805              dj[iSequence]-=difference;
3806              cost[iSequence]-=difference;
3807            }
3808#endif
3809          }
3810          abcMatrix_->add(output, iSequence, movement);
3811          objectiveChange_ += cost[iSequence]*movement;
3812        }
3813      }
3814      array.setNumElements(0);
3815      // double check something flipped
3816      if (output.getNumElements()) {
3817#if ABC_PARALLEL==0
3818        abcFactorization_->updateColumn(output);
3819#else
3820        assert (FACTOR_CPU>2);
3821        abcFactorization_->updateColumnCpu(output,1);
3822#endif
3823        abcDualRowPivot_->updatePrimalSolution(output,1.0);
3824      }
3825#if CLP_CAN_HAVE_ZERO_OBJ>1
3826    } else {
3827      clearArrays(3);
3828    }
3829#endif
3830  }
3831  // recompute dualOut_
3832  valueOut_ = solutionBasic_[pivotRow_];
3833  if (directionOut_ < 0) {
3834    dualOut_ = valueOut_ - upperOut_;
3835  } else {
3836    dualOut_ = lowerOut_ - valueOut_;
3837  }
3838#if 0
3839  // amount primal will move
3840  movement_ = -dualOut_ * directionOut_ / alpha_;
3841  double movementOld = oldDualOut * directionOut_ / alpha_;
3842  // so objective should increase by fabs(dj)*movement_
3843  // but we already have objective change - so check will be good
3844  if (objectiveChange_ + fabs(movementOld * dualIn_) < -CoinMax(1.0e-5, 1.0e-12 * fabs(rawObjectiveValue_))) {
3845#if ABC_NORMAL_DEBUG>3
3846    if (handler_->logLevel() & 32)
3847      printf("movement %g, movementOld %g swap change %g, rest %g  * %g\n",
3848             objectiveChange_ + fabs(movement_ * dualIn_),
3849             movementOld,objectiveChange_, movement_, dualIn_);
3850#endif
3851  }
3852  if (0) {
3853    if(abcFactorization_->pivots()) {
3854      // going backwards - factorize
3855      problemStatus_ = -2; // factorize now
3856      stateOfIteration_=2;
3857    }
3858  }
3859#endif
3860  return numberFlipped;
3861}
3862/* Undo a flip
3863 */
3864void 
3865AbcSimplexDual::flipBack(int number)
3866{
3867  CoinIndexedVector & array = usefulArray_[arrayForFlipBounds_];
3868  const int *  COIN_RESTRICT which=array.getIndices();
3869  double *  COIN_RESTRICT solution = abcSolution_;
3870  const double *  COIN_RESTRICT lower = abcLower_;
3871  const double *  COIN_RESTRICT upper = abcUpper_;
3872  for (int i = 0; i < number; i++) {
3873    int iSequence = which[i];
3874    unsigned char iStatus = internalStatus_[iSequence]&3;
3875    assert ((internalStatus_[iSequence]&7)==0||
3876            (internalStatus_[iSequence]&7)==1);
3877    if (iStatus==0) {
3878      // at lower bound
3879      // to upper bound
3880      setInternalStatus(iSequence, atUpperBound);
3881      solution[iSequence] = upper[iSequence];
3882    } else {
3883      // at upper bound
3884      // to lower bound
3885      setInternalStatus(iSequence, atLowerBound);
3886      solution[iSequence] = lower[iSequence];
3887    }
3888  }
3889}
3890// Restores bound to original bound
3891void
3892AbcSimplexDual::originalBound( int iSequence)
3893{
3894  if (getFakeBound(iSequence) != noFake) {
3895    numberFake_--;;
3896    setFakeBound(iSequence, noFake);
3897    abcLower_[iSequence] = lowerSaved_[iSequence];
3898    abcUpper_[iSequence] = upperSaved_[iSequence];
3899  }
3900}
3901/* As changeBounds but just changes new bounds for a single variable.
3902   Returns true if change */
3903bool
3904AbcSimplexDual::changeBound( int iSequence)
3905{
3906  // old values
3907  double oldLower = abcLower_[iSequence];
3908  double oldUpper = abcUpper_[iSequence];
3909  double value = abcSolution_[iSequence];
3910  bool modified = false;
3911  // original values
3912  double lowerValue = lowerSaved_[iSequence];
3913  double upperValue = upperSaved_[iSequence];
3914  setFakeBound(iSequence,noFake);
3915  if (value == oldLower) {
3916    if (upperValue > oldLower + currentDualBound_) {
3917      abcUpper_[iSequence] = oldLower + currentDualBound_;
3918      setFakeBound(iSequence, upperFake);
3919      assert (abcLower_[iSequence]==lowerSaved_[iSequence]);
3920      assert (abcUpper_[iSequence]<upperSaved_[iSequence]);
3921      modified = true;
3922      numberFake_++;
3923    }
3924  } else if (value == oldUpper) {
3925    if (lowerValue < oldUpper - currentDualBound_) {
3926      abcLower_[iSequence] = oldUpper - currentDualBound_;
3927      setFakeBound(iSequence, lowerFake);
3928      assert (abcLower_[iSequence]>lowerSaved_[iSequence]);
3929      assert (abcUpper_[iSequence]==upperSaved_[iSequence]);
3930      modified = true;
3931      numberFake_++;
3932    }
3933  } else {
3934    assert(value == oldLower || value == oldUpper);
3935  }
3936  return modified;
3937}
3938/* Checks number of variables at fake bounds.  This is used by fastDual
3939   so can exit gracefully before end */
3940int
3941AbcSimplexDual::numberAtFakeBound()
3942{
3943  int numberFake = 0;
3944
3945  if (numberFake_) {
3946    for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
3947      FakeBound bound = getFakeBound(iSequence);
3948      switch(getInternalStatus(iSequence)) {
3949       
3950      case basic:
3951        break;
3952      case isFree:
3953      case superBasic:
3954      case AbcSimplex::isFixed:
3955        //setFakeBound (iSequence, noFake);
3956        break;
3957      case atUpperBound:
3958        if (bound == upperFake || bound == bothFake)
3959          numberFake++;
3960        break;
3961      case atLowerBound:
3962        if (bound == lowerFake || bound == bothFake)
3963          numberFake++;
3964        break;
3965      }
3966    }
3967  }
3968  return numberFake;
3969}
3970int
3971AbcSimplexDual::resetFakeBounds(int type)
3972{
3973  if (type == 0||type==1) {
3974    // put back original bounds and then check
3975    //CoinAbcMemcpy(abcLower_,lowerSaved_,numberTotal_);
3976    //CoinAbcMemcpy(abcUpper_,upperSaved_,numberTotal_);
3977    double dummyChangeCost = 0.0;
3978    return changeBounds(3+type, dummyChangeCost);
3979  } else if (type < 0) {
3980#ifndef NDEBUG
3981    // just check
3982    int iSequence;
3983    int nFake = 0;
3984    int nErrors = 0;
3985    int nSuperBasic = 0;
3986    int nWarnings = 0;
3987    for (iSequence = 0; iSequence < numberTotal_; iSequence++) {
3988      FakeBound fakeStatus = getFakeBound(iSequence);
3989      Status status = getInternalStatus(iSequence);
3990      bool isFake = false;
3991#ifdef CLP_INVESTIGATE
3992      char RC = 'C';
3993      int jSequence = iSequence;
3994      if (jSequence >= numberColumns_) {
3995        RC = 'R';
3996        jSequence -= numberColumns_;
3997      }
3998#endif
3999      double lowerValue = lowerSaved_[iSequence];
4000      double upperValue = upperSaved_[iSequence];
4001      double value = abcSolution_[iSequence];
4002      CoinRelFltEq equal;
4003      if (status == atUpperBound ||
4004          status == atLowerBound) {
4005        if (fakeStatus == AbcSimplexDual::upperFake) {
4006          if(!equal(abcUpper_[iSequence], (lowerValue + currentDualBound_)) ||
4007             !(equal(abcUpper_[iSequence], value) ||
4008               equal(abcLower_[iSequence], value))) {
4009            nErrors++;
4010#ifdef CLP_INVESTIGATE
4011            printf("** upperFake %c%d %g <= %g <= %g true %g, %g\n",
4012                   RC, jSequence, abcLower_[iSequence], abcSolution_[iSequence],
4013                   abcUpper_[iSequence], lowerValue, upperValue);
4014#endif
4015          }
4016          isFake = true;;
4017        } else if (fakeStatus == AbcSimplexDual::lowerFake) {
4018          if(!equal(abcLower_[iSequence], (upperValue - currentDualBound_)) ||
4019             !(equal(abcUpper_[iSequence], value) ||
4020               equal(abcLower_[iSequence], value))) {
4021            nErrors++;
4022#ifdef CLP_INVESTIGATE
4023            printf("** lowerFake %c%d %g <= %g <= %g true %g, %g\n",
4024                   RC, jSequence, abcLower_[iSequence], abcSolution_[iSequence],
4025                   abcUpper_[iSequence], lowerValue, upperValue);
4026#endif
4027          }
4028          isFake = true;;
4029        } else if (fakeStatus == AbcSimplexDual::bothFake) {
4030          nWarnings++;
4031#ifdef CLP_INVESTIGATE
4032          printf("** %d at bothFake?\n", iSequence);
4033#endif
4034        } else if (abcUpper_[iSequence] - abcLower_[iSequence] > 2.0 * currentDualBound_) {
4035          nErrors++;
4036#ifdef CLP_INVESTIGATE
4037          printf("** noFake! %c%d %g <= %g <= %g true %g, %g\n",
4038                 RC, jSequence, abcLower_[iSequence], abcSolution_[iSequence],
4039                 abcUpper_[iSequence], lowerValue, upperValue);
4040#endif
4041        }
4042      } else if (status == superBasic || status == isFree) {
4043        nSuperBasic++;
4044        //printf("** free or superbasic %c%d %g <= %g <= %g true %g, %g - status %d\n",
4045        //     RC,jSequence,abcLower_[iSequence],abcSolution_[iSequence],
4046        //     abcUpper_[iSequence],lowerValue,upperValue,status);
4047      } else if (status == basic) {
4048        bool odd = false;
4049        if (!equal(abcLower_[iSequence], lowerValue))
4050          odd = true;
4051        if (!equal(abcUpper_[iSequence], upperValue))
4052          odd = true;
4053        if (odd) {
4054#ifdef CLP_INVESTIGATE
4055          printf("** basic %c%d %g <= %g <= %g true %g, %g\n",
4056                 RC, jSequence, abcLower_[iSequence], abcSolution_[iSequence],
4057                 abcUpper_[iSequence], lowerValue, upperValue);
4058#endif
4059          nWarnings++;
4060        }
4061      } else if (status == isFixed) {
4062        if (!equal(abcUpper_[iSequence], abcLower_[iSequence])) {
4063          nErrors++;
4064#ifdef CLP_INVESTIGATE
4065          printf("** fixed! %c%d %g <= %g <= %g true %g, %g\n",
4066                 RC, jSequence, abcLower_[iSequence], abcSolution_[iSequence],
4067                 abcUpper_[iSequence], lowerValue, upperValue);
4068#endif
4069        }
4070      }
4071      if (isFake) {
4072        nFake++;
4073      } else {
4074        if (fakeStatus != AbcSimplexDual::noFake) {
4075          nErrors++;
4076#ifdef CLP_INVESTIGATE
4077          printf("** bad fake status %c%d %d\n",
4078                 RC, jSequence, fakeStatus);
4079#endif
4080        }
4081      }
4082    }
4083    if (nFake != numberFake_) {
4084#ifdef CLP_INVESTIGATE
4085      printf("nfake %d numberFake %d\n", nFake, numberFake_);
4086#endif
4087      nErrors++;
4088    }
4089    if (nErrors || type <= -1000) {
4090#ifdef CLP_INVESTIGATE
4091      printf("%d errors, %d warnings, %d free/superbasic, %d fake\n",
4092             nErrors, nWarnings, nSuperBasic, numberFake_);
4093      printf("currentDualBound %g\n",
4094             currentDualBound_);
4095#endif
4096      if (type <= -1000) {
4097        iSequence = -type;
4098        iSequence -= 1000;
4099#ifdef CLP_INVESTIGATE
4100        char RC = 'C';
4101        int jSequence = iSequence;
4102        if (jSequence >= numberColumns_) {
4103          RC = 'R';
4104          jSequence -= numberColumns_;
4105        }
4106        double lowerValue = lowerSaved_[iSequence];
4107        double upperValue = upperSaved_[iSequence];
4108        printf("*** movement>1.0e30 for  %c%d %g <= %g <= %g true %g, %g - status %d\n",
4109               RC, jSequence, abcLower_[iSequence], abcSolution_[iSequence],
4110               abcUpper_[iSequence], lowerValue, upperValue, internalStatus_[iSequence]);
4111#endif
4112        assert (nErrors); // should have been picked up
4113      }
4114      assert (!nErrors);
4115    }
4116#endif
4117  } else {
4118    CoinAbcMemcpy(abcLower_,lowerSaved_,numberTotal_);
4119    CoinAbcMemcpy(abcUpper_,upperSaved_,numberTotal_);
4120    // reset using status
4121    numberFake_ = 0;
4122    for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
4123      FakeBound fakeStatus = getFakeBound(iSequence);
4124      if (fakeStatus != AbcSimplexDual::noFake) {
4125        Status status = getInternalStatus(iSequence);
4126        if (status == basic) {
4127          setFakeBound(iSequence, AbcSimplexDual::noFake);
4128          continue;
4129        }
4130        double lowerValue = abcLower_[iSequence];
4131        double upperValue = abcUpper_[iSequence];
4132        double value = abcSolution_[iSequence];
4133        numberFake_++;
4134        if (fakeStatus == AbcSimplexDual::upperFake) {
4135          abcUpper_[iSequence] = lowerValue + currentDualBound_;
4136          if (status == AbcSimplex::atLowerBound) {
4137            abcSolution_[iSequence] = lowerValue;
4138          } else if (status == AbcSimplex::atUpperBound) {
4139            abcSolution_[iSequence] = abcUpper_[iSequence];
4140          } else {
4141            abort();
4142          }
4143        } else if (fakeStatus == AbcSimplexDual::lowerFake) {
4144          abcLower_[iSequence] = upperValue - currentDualBound_;
4145          if (status == AbcSimplex::atLowerBound) {
4146            abcSolution_[iSequence] = abcLower_[iSequence];
4147          } else if (status == AbcSimplex::atUpperBound) {
4148            abcSolution_[iSequence] = upperValue;
4149          } else {
4150            abort();
4151          }
4152        } else {
4153          assert (fakeStatus == AbcSimplexDual::bothFake);
4154          if (status == AbcSimplex::atLowerBound) {
4155            abcLower_[iSequence] = value;
4156            abcUpper_[iSequence] = value + currentDualBound_;
4157          } else if (status == AbcSimplex::atUpperBound) {
4158            abcUpper_[iSequence] = value;
4159            abcLower_[iSequence] = value - currentDualBound_;
4160          } else if (status == AbcSimplex::isFree ||
4161                     status == AbcSimplex::superBasic) {
4162            abcLower_[iSequence] = value - 0.5 * currentDualBound_;
4163            abcUpper_[iSequence] = value + 0.5 * currentDualBound_;
4164          } else {
4165            abort();
4166          }
4167        }
4168      }
4169    }
4170  }
4171  return -1;
4172}
4173// maybe could JUST check against perturbationArray
4174void
4175AbcSimplex::checkDualSolutionPlusFake()
4176{
4177 
4178  sumDualInfeasibilities_ = 0.0;
4179  numberDualInfeasibilities_ = 0;
4180  sumFakeInfeasibilities_=0.0;
4181  int firstFreePrimal = -1;
4182  int firstFreeDual = -1;
4183  int numberSuperBasicWithDj = 0;
4184  bestPossibleImprovement_ = 0.0;
4185  // we can't really trust infeasibilities if there is dual error
4186  double error = CoinMin(1.0e-2, largestDualError_);
4187  // allow tolerance at least slightly bigger than standard
4188  double relaxedTolerance = currentDualTolerance_ +  error;
4189  // allow bigger tolerance for possible improvement
4190  double possTolerance = 5.0 * relaxedTolerance;
4191  sumOfRelaxedDualInfeasibilities_ = 0.0;
4192  //rawObjectiveValue_ -=sumNonBasicCosts_;
4193  sumNonBasicCosts_ = 0.0;
4194  const double *  COIN_RESTRICT abcLower = abcLower_;
4195  const double *  COIN_RESTRICT abcUpper = abcUpper_;
4196  const double *  COIN_RESTRICT abcSolution = abcSolution_;
4197#ifdef CLEANUP_DJS
4198  double *  COIN_RESTRICT abcDj = abcDj_;
4199  double *  COIN_RESTRICT abcCost = abcCost_;
4200  // see if we can clean up djs
4201  // may have perturbation
4202  const double * COIN_RESTRICT originalCost = abcPerturbation_;
4203#else
4204  const double *  COIN_RESTRICT abcDj = abcDj_;
4205  const double *  COIN_RESTRICT abcCost = abcCost_;
4206#endif
4207  //const double *  COIN_RESTRICT abcPerturbation = abcPerturbation_;
4208  int badDjSequence=-1;
4209  for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
4210    if (getInternalStatus(iSequence) != basic) {
4211      sumNonBasicCosts_ += abcSolution[iSequence]*abcCost[iSequence];
4212      // not basic
4213      double distanceUp = abcUpper[iSequence] - abcSolution[iSequence];
4214      double distanceDown = abcSolution[iSequence] - abcLower[iSequence];
4215      double value = abcDj[iSequence];
4216      if (distanceUp > primalTolerance_) {
4217        // Check if "free"
4218        if (distanceDown <= primalTolerance_) {
4219          // should not be negative
4220          if (value < 0.0) {
4221            value = - value;
4222            if (value > currentDualTolerance_) {
4223              bool treatAsSuperBasic=false;
4224              if (getInternalStatus(iSequence)==superBasic) {
4225#ifdef PAN
4226                if (fakeSuperBasic(iSequence)>=0) {
4227#endif
4228                  treatAsSuperBasic=true;
4229                  numberSuperBasicWithDj++;
4230                  if (firstFreeDual < 0)
4231                    firstFreeDual = iSequence;
4232#ifdef PAN
4233                }
4234#endif
4235              }
4236              if (!treatAsSuperBasic) {
4237                sumDualInfeasibilities_ += value - currentDualTolerance_;
4238                badDjSequence=iSequence;
4239                if (value > possTolerance)
4240                  bestPossibleImprovement_ += CoinMin(distanceUp, 1.0e10) * value;
4241                if (value > relaxedTolerance)
4242                  sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
4243                numberDualInfeasibilities_ ++;
4244              }
4245            }
4246#ifdef CLEANUP_DJS
4247          } else {
4248            // adjust cost etc
4249            if (abcCost[iSequence]>originalCost[iSequence]) {
4250              double difference = CoinMin(abcCost[iSequence]-originalCost[iSequence],value);
4251              //printf("Lb %d changing cost from %g to %g, dj from %g to %g\n",
4252              //     iSequence,abcCost[iSequence],abcCost[iSequence]-difference,
4253              //     abcDj[iSequence],abcDj[iSequence]-difference);
4254              abcDj[iSequence]=value-difference;
4255              abcCost[iSequence]-=difference;
4256            }
4257#endif
4258          }
4259        } else {
4260          // free
4261          value=fabs(value);
4262          if (value > 1.0 * relaxedTolerance) {
4263            numberSuperBasicWithDj++;
4264            if (firstFreeDual < 0)
4265              firstFreeDual = iSequence;
4266          }
4267          if (value > currentDualTolerance_) {
4268            sumDualInfeasibilities_ += value - currentDualTolerance_;
4269            if (value > possTolerance)
4270              bestPossibleImprovement_ += CoinMin(distanceUp, 1.0e10) * value;
4271            if (value > relaxedTolerance)
4272              sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
4273            numberDualInfeasibilities_ ++;
4274          }
4275          if (firstFreePrimal < 0)
4276            firstFreePrimal = iSequence;
4277        }
4278      } else if (distanceDown > primalTolerance_) {
4279        // should not be positive
4280        if (value > 0.0) {
4281          if (value > currentDualTolerance_) {
4282            bool treatAsSuperBasic=false;
4283            if (getInternalStatus(iSequence)==superBasic) {
4284#ifdef PAN
4285              if (fakeSuperBasic(iSequence)>=0) {
4286#endif
4287                treatAsSuperBasic=true;
4288                numberSuperBasicWithDj++;
4289                if (firstFreeDual < 0)
4290                  firstFreeDual = iSequence;
4291#ifdef PAN
4292              }
4293#endif
4294            }
4295            if (!treatAsSuperBasic) {
4296              sumDualInfeasibilities_ += value - currentDualTolerance_;
4297              badDjSequence=iSequence;
4298              if (value > possTolerance)
4299                bestPossibleImprovement_ += value * CoinMin(distanceDown, 1.0e10);
4300              if (value > relaxedTolerance)
4301                sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
4302              numberDualInfeasibilities_ ++;
4303            }
4304          }
4305#ifdef CLEANUP_DJS
4306        } else {
4307          // adjust cost etc
4308          if (abcCost[iSequence]<originalCost[iSequence]) {
4309            double difference = CoinMax(abcCost[iSequence]-originalCost[iSequence],value);
4310            //printf("Ub %d changing cost from %g to %g, dj from %g to %g\n",
4311            //       iSequence,abcCost[iSequence],abcCost[iSequence]-difference,
4312            //       abcDj[iSequence],abcDj[iSequence]-difference);
4313            abcDj[iSequence]=value-difference;
4314            abcCost[iSequence]-=difference;
4315          }
4316#endif
4317        }
4318      }
4319    }
4320  }
4321#ifdef CLEANUP_DJS
4322  if (algorithm_<0&&numberDualInfeasibilities_==1&&firstFreeDual<0&&sumDualInfeasibilities_<1.0e-1) {
4323    // clean bad one
4324    assert (badDjSequence>=0);
4325    abcCost[badDjSequence]-=abcDj_[badDjSequence];
4326    abcDj_[badDjSequence]=0.0;
4327  }
4328#endif
4329  numberDualInfeasibilitiesWithoutFree_ = numberDualInfeasibilities_-numberSuperBasicWithDj;
4330  numberFreeNonBasic_=numberSuperBasicWithDj;
4331  firstFree_=-1;
4332  if (algorithm_ < 0 && firstFreeDual >= 0) {
4333    // dual
4334    firstFree_ = firstFreeDual;
4335  } else if ((numberSuperBasicWithDj ||
4336              (abcProgress_.lastIterationNumber(0) <= 0))&&algorithm_>0) {
4337    firstFree_ = firstFreePrimal;
4338  }
4339}
4340// Perturbs problem
4341void 
4342AbcSimplexDual::perturb(double factor)
4343{
4344  const double * COIN_RESTRICT lower = lowerSaved_;
4345  const double * COIN_RESTRICT upper = upperSaved_;
4346  const double * COIN_RESTRICT perturbation = abcPerturbation_;
4347  double * COIN_RESTRICT cost = abcCost_;
4348  for (int i=0;i<numberTotal_;i++) {
4349    if (getInternalStatus(i)!=basic) {
4350      if (fabs(lower[i])<fabs(upper[i])) {
4351        cost[i]+=factor*perturbation[i];
4352      } else if (upper[i]!=COIN_DBL_MAX&&upper[i]!=lower[i]) {
4353        cost[i]-=factor*perturbation[i];
4354      } else {
4355        // free
4356        //cost[i]=costSaved[i];
4357      }
4358    }
4359  }
4360#ifdef CLEANUP_DJS
4361  // copy cost to perturbation
4362  CoinAbcMemcpy(abcPerturbation_,abcCost_,numberTotal_);
4363#endif
4364}
4365#if ABC_NORMAL_DEBUG>0
4366#define PERT_STATISTICS
4367#endif
4368#ifdef PERT_STATISTICS
4369static void breakdown(const char * name, int numberLook, const double * region)
4370{
4371     double range[] = {
4372          -COIN_DBL_MAX,
4373          -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,
4374          -1.0, -3.0e-1,
4375          -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,
4376          0.0,
4377          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,
4378          1.0, 3.0,
4379          1.0e1, 3.0e1, 1.0e2, 3.0e2, 1.0e3, 3.0e3, 1.0e4, 3.0e4, 1.0e5, 1.0e8, 1.0e11, 1.0e15,
4380          COIN_DBL_MAX
4381     };
4382     int nRanges = static_cast<int> (sizeof(range) / sizeof(double));
4383     int * number = new int[nRanges];
4384     memset(number, 0, nRanges * sizeof(int));
4385     int * numberExact = new int[nRanges];
4386     memset(numberExact, 0, nRanges * sizeof(int));
4387     int i;
4388     for ( i = 0; i < numberLook; i++) {
4389          double value = region[i];
4390          for (int j = 0; j < nRanges; j++) {
4391               if (value == range[j]) {
4392                    numberExact[j]++;
4393                    break;
4394               } else if (value < range[j]) {
4395                    number[j]++;
4396                    break;
4397               }
4398          }
4399     }
4400     printf("\n%s has %d entries\n", name, numberLook);
4401     for (i = 0; i < nRanges; i++) {
4402          if (number[i])
4403               printf("%d between %g and %g", number[i], range[i-1], range[i]);
4404          if (numberExact[i]) {
4405               if (number[i])
4406                    printf(", ");
4407               printf("%d exactly at %g", numberExact[i], range[i]);
4408          }
4409          if (number[i] + numberExact[i])
4410               printf("\n");
4411     }
4412     delete [] number;
4413     delete [] numberExact;
4414}
4415#endif
4416// Perturbs problem
4417void 
4418AbcSimplexDual::perturbB(double /*factorIn*/,int /*type*/)
4419{
4420  double overallMultiplier= (perturbation_==50||perturbation_>54) ? 1.0 : 0.1;
4421  // dual perturbation
4422  double perturbation = 1.0e-20;
4423  // maximum fraction of cost to perturb
4424  double maximumFraction = 1.0e-5;
4425  int maxLength = 0;
4426  int minLength = numberRows_;
4427  double averageCost = 0.0;
4428  int numberNonZero = 0;
4429  // See if we need to perturb
4430  double * COIN_RESTRICT sort = new double[numberColumns_];
4431  // Use objective BEFORE scaling ??
4432  const double * COIN_RESTRICT obj = objective();
4433  for (int i = 0; i < numberColumns_; i++) {
4434    double value = fabs(obj[i]);
4435    sort[i] = value;
4436    averageCost += fabs(abcCost_[i+maximumAbcNumberRows_]);
4437    if (value)
4438      numberNonZero++;
4439  }
4440  if (numberNonZero)
4441    averageCost /= static_cast<double> (numberNonZero);
4442  else
4443    averageCost = 1.0;
4444  // allow for cost scaling
4445  averageCost *= objectiveScale_;
4446  std::sort(sort, sort + numberColumns_);
4447  int number = 1;
4448  double last = sort[0];
4449  for (int i = 1; i < numberColumns_; i++) {
4450    if (last != sort[i])
4451      number++;
4452    last = sort[i];
4453  }
4454  delete [] sort;
4455#ifdef PERT_STATISTICS
4456  printf("%d different values\n",number);
4457#endif
4458#if 0
4459  printf("nnz %d percent %d", number, (number * 100) / numberColumns_);
4460  if (number * 4 > numberColumns_)
4461    printf(" - Would not perturb\n");
4462  else
4463    printf(" - Would perturb\n");
4464#endif
4465  // If small costs then more dangerous (better to scale)
4466  double numberMultiplier=1.25;
4467  if (averageCost<1.0e-3)
4468    numberMultiplier=20.0;
4469  else if (averageCost<1.0e-2)
4470    numberMultiplier=10.0;
4471  else if (averageCost<1.0e-1)
4472    numberMultiplier=3.0;
4473  else if (averageCost<1.0)
4474    numberMultiplier=2.0;
4475  numberMultiplier *= 2.0; // try
4476  if (number * numberMultiplier > numberColumns_||perturbation_>=100) {
4477    perturbation_=101;
4478#if ABC_NORMAL_DEBUG>0
4479    printf("Not perturbing - %d different costs, average %g\n",
4480           number,averageCost);
4481#endif
4482#ifdef CLEANUP_DJS
4483    // copy cost to perturbation
4484    CoinAbcMemcpy(abcPerturbation_,abcCost_,numberTotal_);
4485#endif
4486    return ; // good enough
4487  }
4488  const int * COIN_RESTRICT columnLength = matrix_->getVectorLengths()-maximumAbcNumberRows_;
4489  const double * COIN_RESTRICT lower = lowerSaved_;
4490  const double * COIN_RESTRICT upper = upperSaved_;
4491  for (int iSequence = maximumAbcNumberRows_; iSequence < numberTotal_; iSequence++) {
4492    if (lower[iSequence] < upper[iSequence]) {
4493      int length = columnLength[iSequence];
4494      if (length > 2) {
4495        maxLength = CoinMax(maxLength, length);
4496        minLength = CoinMin(minLength, length);
4497      }
4498    }
4499  }
4500  bool uniformChange = false;
4501  bool inCbcOrOther = (specialOptions_ & 0x03000000) != 0;
4502  double constantPerturbation = 10.0 * dualTolerance_;
4503#ifdef PERT_STATISTICS
4504  breakdown("Objective before perturbation", numberColumns_, abcCost_+numberRows_);
4505#endif
4506  //#define PERTURB_OLD_WAY
4507#ifndef PERTURB_OLD_WAY
4508#if 1
4509#ifdef HEAVY_PERTURBATION
4510    if (perturbation_==50)
4511      perturbation_=HEAVY_PERTURBATION;
4512#endif
4513  if (perturbation_ >= 50) {
4514    // Experiment
4515    // maximumFraction could be 1.0e-10 to 1.0
4516    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};
4517    int whichOne = CoinMin(perturbation_ - 51,10);
4518    if (whichOne<0) {
4519      if (numberRows_<5000)
4520        whichOne=2; // treat 50 as if 53
4521      else
4522        whichOne=2; // treat 50 as if 52
4523    }
4524    if (inCbcOrOther&&whichOne>0)
4525      whichOne--;
4526    maximumFraction = m[whichOne];
4527    if (whichOne>7)
4528      constantPerturbation *= 10.0; 
4529    //if (whichOne<2)
4530    //constantPerturbation = maximumFraction*1.0e-10 * dualTolerance_;
4531  } else if (inCbcOrOther) {
4532    maximumFraction = 1.0e-6;
4533  }
4534#endif
4535  double smallestNonZero = 1.0e100;
4536  bool modifyRowCosts=false;
4537  numberNonZero = 0;
4538  //perturbation = 1.0e-10;
4539  perturbation = CoinMax(1.0e-10,maximumFraction);
4540  double * COIN_RESTRICT cost = abcCost_;
4541  bool allSame = true;
4542  double lastValue = 0.0;
4543 
4544  for (int iRow = 0; iRow < numberRows_; iRow++) {
4545    double lo = lower[iRow];
4546    double up = upper[iRow];
4547    if (lo < up) {
4548      double value = fabs(cost[iRow]);
4549      perturbation = CoinMax(perturbation, value);
4550      if (value) {
4551        modifyRowCosts = true;
4552        smallestNonZero = CoinMin(smallestNonZero, value);
4553      }
4554    }
4555    if (lo && lo > -1.0e10) {
4556      numberNonZero++;
4557      lo = fabs(lo);
4558      if (!lastValue)
4559        lastValue = lo;
4560      else if (fabs(lo - lastValue) > 1.0e-7)
4561        allSame = false;
4562    }
4563    if (up && up < 1.0e10) {
4564      numberNonZero++;
4565      up = fabs(up);
4566      if (!lastValue)
4567        lastValue = up;
4568      else if (fabs(up - lastValue) > 1.0e-7)
4569        allSame = false;
4570    }
4571  }
4572  double lastValue2 = 0.0;
4573  for (int iSequence = maximumAbcNumberRows_; iSequence < numberTotal_; iSequence++) {
4574    double lo = lower[iSequence];
4575    double up = upper[iSequence];
4576    if (lo < up) {
4577      double value =
4578        fabs(cost[iSequence]);
4579      perturbation = CoinMax(perturbation, value);
4580      if (value) {
4581        smallestNonZero = CoinMin(smallestNonZero, value);
4582      }
4583    }
4584    if (lo && lo > -1.0e10) {
4585      //numberNonZero++;
4586      lo = fabs(lo);
4587      if (!lastValue2)
4588        lastValue2 = lo;
4589      else if (fabs(lo - lastValue2) > 1.0e-7)
4590        allSame = false;
4591    }
4592    if (up && up < 1.0e10) {
4593      //numberNonZero++;
4594      up = fabs(up);
4595      if (!lastValue2)
4596        lastValue2 = up;
4597      else if (fabs(up - lastValue2) > 1.0e-7)
4598        allSame = false;
4599    }
4600  }
4601  if (allSame) {
4602    // Check elements
4603    double smallestNegative;
4604    double largestNegative;
4605    double smallestPositive;
4606    double largestPositive;
4607    matrix_->rangeOfElements(smallestNegative, largestNegative,
4608                             smallestPositive, largestPositive);
4609    if (smallestNegative == largestNegative &&
4610        smallestPositive == largestPositive) {
4611      // Really hit perturbation
4612      double adjust = CoinMin(100.0 * maximumFraction, 1.0e-3 * CoinMax(lastValue, lastValue2));
4613      maximumFraction = CoinMax(adjust, maximumFraction);
4614    }
4615  }
4616  perturbation = CoinMin(perturbation, smallestNonZero / maximumFraction);
4617  double largestZero = 0.0;
4618  double largest = 0.0;
4619  double largestPerCent = 0.0;
4620  // modify costs
4621  bool printOut = (handler_->logLevel() == 63);
4622  printOut = false;
4623  //assert (!modifyRowCosts);
4624  modifyRowCosts = false;
4625  const double * COIN_RESTRICT perturbationArray = perturbationSaved_;
4626  if (modifyRowCosts) {
4627    for (int iRow = 0; iRow < numberRows_; iRow++) {
4628      if (lower[iRow] < upper[iRow]) {
4629        double value = perturbation;
4630        double currentValue = cost[iRow];
4631        value = CoinMin(value, maximumFraction * (fabs(currentValue) + 1.0e-1 * perturbation + 1.0e-3));
4632        double perturbationValue=2.0*(perturbationArray[iRow]-0.5); // for now back to normal random
4633        if (lower[iRow] > -largeValue_) {
4634          if (fabs(lower[iRow]) < fabs(upper[iRow]))
4635            value *= perturbationValue;
4636          else
4637            value *= -perturbationValue;
4638        } else if (upper[iRow] < largeValue_) {
4639          value *= -perturbationValue;
4640        } else {
4641          value = 0.0;
4642        }
4643        if (currentValue) {
4644          largest = CoinMax(largest, fabs(value));
4645          if (fabs(value) > fabs(currentValue)*largestPerCent)
4646            largestPerCent = fabs(value / currentValue);
4647        } else {
4648          largestZero = CoinMax(largestZero, fabs(value));
4649        }
4650        if (printOut)
4651          printf("row %d cost %g change %g\n", iRow, cost[iRow], value);
4652        cost[iRow] += value;
4653      }
4654    }
4655  }
4656
4657  //double extraWeight = 10.0;
4658  // 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};
4659  // 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};
4660  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};
4661  //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};
4662  // Scale back if wanted
4663  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};
4664  if (constantPerturbation < 99.0 * dualTolerance_&&false) {
4665    perturbation *= 0.1;
4666    //extraWeight = 0.5;
4667    memcpy(weight, weight2, sizeof(weight2));
4668  }
4669  // adjust weights if all columns long
4670  double factor = 1.0;
4671  if (maxLength) {
4672    factor = 3.0 / static_cast<double> (minLength);
4673  }
4674  // Make variables with more elements more expensive
4675  const double m1 = 0.5;
4676  double smallestAllowed = overallMultiplier*CoinMin(1.0e-2 * dualTolerance_, maximumFraction);
4677  double largestAllowed = overallMultiplier*CoinMax(1.0e3 * dualTolerance_, maximumFraction * averageCost);
4678  // smaller if in BAB
4679  //if (inCbcOrOther)
4680  //largestAllowed=CoinMin(largestAllowed,1.0e-5);
4681  //smallestAllowed = CoinMin(smallestAllowed,0.1*largestAllowed);
4682  randomNumberGenerator_.setSeed(1000000000);
4683#if ABC_NORMAL_DEBUG>0
4684  printf("perturbation %g constantPerturbation %g smallestNonZero %g maximumFraction %g smallestAllowed %g\n",
4685         perturbation,constantPerturbation,smallestNonZero,maximumFraction,
4686         smallestAllowed);
4687#endif
4688  for (int iSequence = maximumAbcNumberRows_; iSequence < numberTotal_; iSequence++) {
4689    if (lower[iSequence] < upper[iSequence] && getInternalStatus(iSequence) != basic) {
4690      double value = perturbation;
4691      double currentValue = cost[iSequence];
4692      double currentLargestAllowed=fabs(currentValue)*1.0e-3;
4693      if (!currentValue)
4694        currentLargestAllowed=1.0e-3;
4695      if (currentLargestAllowed<=10.0*smallestAllowed)
4696        continue; // don't perturb tiny value
4697      value = CoinMin(value, constantPerturbation + maximumFraction * (fabs(currentValue) + 1.0e-1 * perturbation + 1.0e-8));
4698      //value = CoinMin(value,constantPerturbation;+maximumFraction*fabs(currentValue));
4699      double value2 = constantPerturbation + 1.0e-1 * smallestNonZero;
4700      if (uniformChange) {
4701        value = maximumFraction;
4702        value2 = maximumFraction;
4703      }
4704      //double perturbationValue=2.0*(perturbationArray[iSequence]-0.5); // for now back to normal random
4705      double perturbationValue=randomNumberGenerator_.randomDouble();
4706      perturbationValue = 1.0-m1+m1*perturbationValue; // clean up
4707      if (lower[iSequence] > -largeValue_) {
4708        if (fabs(lower[iSequence]) <
4709            fabs(upper[iSequence])) {
4710          value *= perturbationValue;
4711          value2 *= perturbationValue;
4712        } else {
4713          //value *= -(1.0-m1+m1*randomNumberGenerator_.randomDouble());
4714          //value2 *= -(1.0-m1+m1*randomNumberGenerator_.randomDouble());
4715          value = 0.0;
4716        }
4717      } else if (upper[iSequence] < largeValue_) {
4718        value *= -perturbationValue;
4719        value2 *= -perturbationValue;
4720      } else {
4721        value = 0.0;
4722      }
4723      if (value) {
4724        int length = columnLength[iSequence];
4725        if (length > 3) {
4726          length = static_cast<int> (static_cast<double> (length) * factor);
4727          length = CoinMax(3, length);
4728        }
4729        double multiplier;
4730        if (length < 10)
4731          multiplier = weight[length];
4732        else
4733          multiplier = weight[10];
4734        value *= multiplier;
4735        value = CoinMin(value, value2);
4736        if (value) {
4737          // get in range
4738          if (fabs(value) <= smallestAllowed) {
4739            value *= 10.0;
4740            while (fabs(value) <= smallestAllowed)
4741              value *= 10.0;
4742          } else if (fabs(value) > currentLargestAllowed) {
4743            value *= 0.1;
4744            while (fabs(value) > currentLargestAllowed)
4745              value *= 0.1;
4746          }
4747        }
4748        if (currentValue) {
4749          largest = CoinMax(largest, fabs(value));
4750          if (fabs(value) > fabs(currentValue)*largestPerCent)
4751            largestPerCent = fabs(value / currentValue);
4752        } else {
4753          largestZero = CoinMax(largestZero, fabs(value));
4754        }
4755        // but negative if at ub
4756        if (getInternalStatus(iSequence) == atUpperBound)
4757          value = -value;
4758        if (printOut)
4759          printf("col %d cost %g change %g\n", iSequence, cost[iSequence], value);
4760        cost[iSequence] += value;
4761      }
4762    }
4763  }
4764#else
4765  // old way
4766     if (perturbation_ > 50) {
4767          // Experiment
4768          // maximumFraction could be 1.0e-10 to 1.0
4769          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};
4770          int whichOne = perturbation_ - 51;
4771          //if (inCbcOrOther&&whichOne>0)
4772          //whichOne--;
4773          maximumFraction = m[CoinMin(whichOne, 10)];
4774     } else if (inCbcOrOther) {
4775          //maximumFraction = 1.0e-6;
4776     }
4777     int iRow;
4778     double smallestNonZero = 1.0e100;
4779     numberNonZero = 0;
4780     if (perturbation_ >= 50) {
4781          perturbation = 1.0e-8;
4782          if (perturbation_ > 50 && perturbation_ < 60) 
4783            perturbation = CoinMax(1.0e-8,maximumFraction);
4784          bool allSame = true;
4785          double lastValue = 0.0;
4786          for (iRow = 0; iRow < numberRows_; iRow++) {
4787               double lo = abcLower_[iRow];
4788               double up = abcUpper_[iRow];
4789               if (lo && lo > -1.0e10) {
4790                    numberNonZero++;
4791                    lo = fabs(lo);
4792                    if (!lastValue)
4793                         lastValue = lo;
4794                    else if (fabs(lo - lastValue) > 1.0e-7)
4795                         allSame = false;
4796               }
4797               if (up && up < 1.0e10) {
4798                    numberNonZero++;
4799                    up = fabs(up);
4800                    if (!lastValue)
4801                         lastValue = up;
4802                    else if (fabs(up - lastValue) > 1.0e-7)
4803                         allSame = false;
4804               }
4805          }
4806          double lastValue2 = 0.0;
4807          for (int iColumn = maximumAbcNumberRows_; iColumn < numberTotal_; iColumn++) {
4808               double lo = abcLower_[iColumn];
4809               double up = abcUpper_[iColumn];
4810               if (lo < up) {
4811                    double value =
4812                         fabs(abcCost_[iColumn]);
4813                    perturbation = CoinMax(perturbation, value);
4814                    if (value) {
4815                         smallestNonZero = CoinMin(smallestNonZero, value);
4816                    }
4817               }
4818               if (lo && lo > -1.0e10) {
4819                    //numberNonZero++;
4820                    lo = fabs(lo);
4821                    if (!lastValue2)
4822                         lastValue2 = lo;
4823                    else if (fabs(lo - lastValue2) > 1.0e-7)
4824                         allSame = false;
4825               }
4826               if (up && up < 1.0e10) {
4827                    //numberNonZero++;
4828                    up = fabs(up);
4829                    if (!lastValue2)
4830                         lastValue2 = up;
4831                    else if (fabs(up - lastValue2) > 1.0e-7)
4832                         allSame = false;
4833               }
4834          }
4835          if (allSame) {
4836               // Check elements
4837               double smallestNegative;
4838               double largestNegative;
4839               double smallestPositive;
4840               double largestPositive;
4841               matrix_->rangeOfElements(smallestNegative, largestNegative,
4842                                        smallestPositive, largestPositive);
4843               if (smallestNegative == largestNegative &&
4844                         smallestPositive == largestPositive) {
4845                    // Really hit perturbation
4846                    double adjust = CoinMin(100.0 * maximumFraction, 1.0e-3 * CoinMax(lastValue, lastValue2));
4847                    maximumFraction = CoinMax(adjust, maximumFraction);
4848               }
4849          }
4850          perturbation = CoinMin(perturbation, smallestNonZero / maximumFraction);
4851     }
4852     double largestZero = 0.0;
4853     double largest = 0.0;
4854     double largestPerCent = 0.0;
4855     // modify costs
4856     bool printOut = (handler_->logLevel() == 63);
4857     printOut = false;
4858     // 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};
4859     // 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};
4860     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};
4861     //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};
4862     //double extraWeight = 10.0;
4863     // Scale back if wanted
4864     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};
4865     if (constantPerturbation < 99.0 * dualTolerance_) {
4866          perturbation *= 0.1;
4867          //extraWeight = 0.5;
4868          memcpy(weight, weight2, sizeof(weight2));
4869     }
4870     // adjust weights if all columns long
4871     double factor = 1.0;
4872     if (maxLength) {
4873          factor = 3.0 / static_cast<double> (minLength);
4874     }
4875     // Make variables with more elements more expensive
4876     const double m1 = 0.5;
4877     double smallestAllowed = CoinMin(1.0e-2 * dualTolerance_, maximumFraction);
4878     double largestAllowed = CoinMax(1.0e3 * dualTolerance_, maximumFraction * averageCost);
4879     // smaller if in BAB
4880     //if (inCbcOrOther)
4881     //largestAllowed=CoinMin(largestAllowed,1.0e-5);
4882     //smallestAllowed = CoinMin(smallestAllowed,0.1*largestAllowed);
4883     for (int iColumn = maximumAbcNumberRows_; iColumn < numberTotal_; iColumn++) {
4884          if (abcLower_[iColumn] < abcUpper_[iColumn] 
4885              && getInternalStatus(iColumn) != basic) {
4886               double value = perturbation;
4887               double currentValue = abcCost_[iColumn];
4888               value = CoinMin(value, constantPerturbation + maximumFraction * (fabs(currentValue) + 1.0e-1 * perturbation + 1.0e-8));
4889               //value = CoinMin(value,constantPerturbation;+maximumFraction*fabs(currentValue));
4890               double value2 = constantPerturbation + 1.0e-1 * smallestNonZero;
4891               if (uniformChange) {
4892                    value = maximumFraction;
4893                    value2 = maximumFraction;
4894               }
4895               if (abcLower_[iColumn] > -largeValue_) {
4896                    if (fabs(abcLower_[iColumn]) <
4897                              fabs(abcUpper_[iColumn])) {
4898                         value *= (1.0 - m1 + m1 * randomNumberGenerator_.randomDouble());
4899                         value2 *= (1.0 - m1 + m1 * randomNumberGenerator_.randomDouble());
4900                    } else {
4901                         //value *= -(1.0-m1+m1*randomNumberGenerator_.randomDouble());
4902                         //value2 *= -(1.0-m1+m1*randomNumberGenerator_.randomDouble());
4903                         value = 0.0;
4904                    }
4905               } else if (abcUpper_[iColumn] < largeValue_) {
4906                    value *= -(1.0 - m1 + m1 * randomNumberGenerator_.randomDouble());
4907                    value2 *= -(1.0 - m1 + m1 * randomNumberGenerator_.randomDouble());
4908               } else {
4909                    value = 0.0;
4910               }
4911               if (value) {
4912                    int length = columnLength[iColumn];
4913                    if (length > 3) {
4914                         length = static_cast<int> (static_cast<double> (length) * factor);
4915                         length = CoinMax(3, length);
4916                    }
4917                    double multiplier;
4918#if 1
4919                    if (length < 10)
4920                         multiplier = weight[length];
4921                    else
4922                         multiplier = weight[10];
4923#else
4924                    if (length < 10)
4925                         multiplier = weight[length];
4926                    else
4927                         multiplier = weight[10] + extraWeight * (length - 10);
4928                    multiplier *= 0.5;
4929#endif
4930                    value *= multiplier;
4931                    value = CoinMin(value, value2);
4932                    if (value) {
4933                         // get in range
4934                         if (fabs(value) <= smallestAllowed) {
4935                              value *= 10.0;
4936                              while (fabs(value) <= smallestAllowed)
4937                                   value *= 10.0;
4938                         } else if (fabs(value) > largestAllowed) {
4939                              value *= 0.1;
4940                              while (fabs(value) > largestAllowed)
4941                                   value *= 0.1;
4942                         }
4943                    }
4944                    if (currentValue) {
4945                         largest = CoinMax(largest, fabs(value));
4946                         if (fabs(value) > fabs(currentValue)*largestPerCent)
4947                              largestPerCent = fabs(value / currentValue);
4948                    } else {
4949                         largestZero = CoinMax(largestZero, fabs(value));
4950                    }
4951                    // but negative if at ub
4952                    if (getInternalStatus(iColumn) == atUpperBound)
4953                         value = -value;
4954                    if (printOut)
4955                         printf("col %d cost %g change %g\n", iColumn, objectiveWork_[iColumn], value);
4956                    abcCost_[iColumn] += value;
4957               }
4958          }
4959     }
4960#endif
4961#ifdef PERT_STATISTICS
4962  {
4963    double averageCost = 0.0;
4964    int numberNonZero = 0;
4965    double * COIN_RESTRICT sort = new double[numberColumns_];
4966    for (int i = 0; i < numberColumns_; i++) {
4967      double value = fabs(abcCost_[i+maximumAbcNumberRows_]);
4968      sort[i] = value;
4969      averageCost += value;
4970      if (value)
4971        numberNonZero++;
4972    }
4973    if (numberNonZero)
4974      averageCost /= static_cast<double> (numberNonZero);
4975    else
4976      averageCost = 1.0;
4977    std::sort(sort, sort + numberColumns_);
4978    int number = 1;
4979    double last = sort[0];
4980    for (int i = 1; i < numberColumns_; i++) {
4981      if (last != sort[i])
4982        number++;
4983    last = sort[i];
4984    }
4985    delete [] sort;
4986#if ABC_NORMAL_DEBUG>0
4987    printf("nnz %d percent %d", number, (number * 100) / numberColumns_);
4988    breakdown("Objective after perturbation", numberColumns_, abcCost_+numberRows_);
4989#endif
4990  }
4991#endif
4992  perturbation_=100;
4993  // copy cost to  perturbation
4994#ifdef CLEANUP_DJS
4995  CoinAbcMemcpy(abcPerturbation_,abcCost_,numberTotal_);
4996#endif
4997  handler_->message(CLP_SIMPLEX_PERTURB, messages_)
4998    << 100.0 * maximumFraction << perturbation << largest << 100.0 * largestPerCent << largestZero
4999    << CoinMessageEol;
5000}
5001/* Does something about fake tolerances
5002   -1 initializes
5003   1 might be optimal (back to original costs)
5004   2 after change currentDualBound
5005   3 might be optimal BUT change currentDualBound
5006*/
5007int 
5008AbcSimplexDual::bounceTolerances(int type)
5009{
5010  // for now turn off
5011  if (type==-1) {
5012    if (!perturbationFactor_)
5013      perturbationFactor_=0.001;
5014    CoinAbcMemcpy(abcCost_,costSaved_,numberTotal_);
5015    if (!perturbationFactor_) {
5016      //currentDualBound_=CoinMin(dualBound_,1.0e11);
5017      //CoinAbcMultiplyAdd(perturbationSaved_,numberTotal_,currentDualTolerance_,
5018      //                 abcPerturbation_,0.0);
5019    } else {
5020      //currentDualBound_=CoinMin(dualBound_,1.0e11);
5021      double factor=currentDualTolerance_*perturbationFactor_;
5022      perturbB(factor,0);
5023    }
5024#if ABC_NORMAL_DEBUG>0
5025    //printf("new dual bound %g (bounceTolerances)\n",currentDualBound_);
5026#endif
5027  } else {
5028    if (stateDualColumn_==-2) {
5029      double newTolerance=CoinMin(1.0e-1,dualTolerance_*100.0);
5030#if ABC_NORMAL_DEBUG>0
5031      printf("After %d iterations reducing dual tolerance from %g to %g\n",
5032             numberIterations_,currentDualTolerance_,newTolerance);
5033#endif
5034      currentDualTolerance_=newTolerance;
5035      stateDualColumn_=-1;
5036    } else if (stateDualColumn_==-1) {
5037#if ABC_NORMAL_DEBUG>0
5038      printf("changing dual tolerance from %g to %g (vaguely optimal)\n",
5039             currentDualTolerance_,dualTolerance_);
5040#endif
5041      currentDualTolerance_=dualTolerance_;
5042      stateDualColumn_=0;
5043    } else {
5044      if (stateDualColumn_>1) {
5045        normalDualColumnIteration_ = numberIterations_+CoinMax(10,numberRows_/4);
5046#if ABC_NORMAL_DEBUG>0
5047        printf("no cost modification until iteration %d\n",normalDualColumnIteration_);
5048#endif
5049      }
5050    }
5051    resetFakeBounds(0);
5052    CoinAbcMemcpy(abcCost_,abcPerturbation_,numberTotal_);
5053    numberChanged_ = 0; // Number of variables with changed costs
5054    moveToBasic();
5055    if (type!=1)
5056      return 0;
5057    // save status and solution
5058    CoinAbcMemcpy(status_,internalStatus_,numberTotal_);
5059    CoinAbcMemcpy(solutionSaved_,abcSolution_,numberTotal_);
5060#if ABC_NORMAL_DEBUG>0
5061    printf("A sumDual %g sumPrimal %g\n",sumDualInfeasibilities_,sumPrimalInfeasibilities_);
5062#endif
5063    gutsOfSolution(3);
5064#if ABC_NORMAL_DEBUG>0
5065    printf("B sumDual %g sumPrimal %g\n",sumDualInfeasibilities_,sumPrimalInfeasibilities_);
5066#endif
5067    double saveSum1=sumPrimalInfeasibilities_;
5068    double saveSumD1=sumDualInfeasibilities_;
5069    bool goToPrimal=false;
5070    if (!sumOfRelaxedDualInfeasibilities_) {
5071      if (perturbation_<101) {
5072        // try going back to original costs
5073        CoinAbcMemcpy(abcCost_,costSaved_,numberTotal_);
5074        perturbation_=101;
5075        numberChanged_ = 0; // Number of variables with changed costs
5076        moveToBasic();
5077        gutsOfSolution(3);
5078#if ABC_NORMAL_DEBUG>0
5079        printf("B2 sumDual %g sumPrimal %g\n",sumDualInfeasibilities_,sumPrimalInfeasibilities_);
5080#endif
5081        saveSumD1=sumDualInfeasibilities_;
5082      }
5083      if (sumOfRelaxedDualInfeasibilities_) {
5084        makeNonFreeVariablesDualFeasible();
5085        moveToBasic();
5086        abcProgress_.reset();
5087        gutsOfSolution(3);
5088#if ABC_NORMAL_DEBUG>0
5089        printf("C sumDual %g sumPrimal %g\n",sumDualInfeasibilities_,sumPrimalInfeasibilities_);
5090#endif
5091        if (sumPrimalInfeasibilities_>1.0e2*saveSum1+currentDualBound_*0.2
5092          /*|| sumOfRelaxedDualInfeasibilities_<1.0e-4*sumPrimalInfeasibilities_*/) {
5093          numberDisasters_++;
5094          if (numberDisasters_>2)
5095            goToPrimal=true;
5096        }
5097      } else {
5098        numberDualInfeasibilities_=0; 
5099        sumDualInfeasibilities_=0.0; 
5100      }
5101    } else {
5102      // carry on
5103      makeNonFreeVariablesDualFeasible();
5104      moveToBasic();
5105      abcProgress_.reset();
5106      gutsOfSolution(3);
5107#if ABC_NORMAL_DEBUG>0
5108      printf("C2 sumDual %g sumPrimal %g\n",sumDualInfeasibilities_,sumPrimalInfeasibilities_);
5109#endif
5110      numberDisasters_++;
5111      if (sumPrimalInfeasibilities_>1.0e2*saveSum1+currentDualBound_*0.2
5112          ||(saveSumD1<1.0e-4*sumPrimalInfeasibilities_&&saveSumD1<1.0)) {
5113        numberDisasters_++;
5114        if (numberDisasters_>2||saveSumD1<1.0e-10*sumPrimalInfeasibilities_)
5115          goToPrimal=true;
5116      }
5117    }
5118    if (goToPrimal) {
5119      // go to primal
5120      CoinAbcMemcpy(internalStatus_,status_,numberTotal_);
5121      CoinAbcMemcpy(abcSolution_,solutionSaved_,numberTotal_);
5122      int *  COIN_RESTRICT abcPivotVariable = abcPivotVariable_;
5123      // redo pivotVariable
5124      int numberBasic=0;
5125      for (int i=0;i<numberTotal_;i++) {
5126        if (getInternalStatus(i)==basic)
5127          abcPivotVariable[numberBasic++]=i;
5128      }
5129      assert (numberBasic==numberRows_);
5130      moveToBasic();
5131      checkPrimalSolution(false);
5132#if ABC_NORMAL_DEBUG>0
5133      printf("D sumPrimal %g\n",sumPrimalInfeasibilities_);
5134#endif
5135      return 1;
5136    }
5137  }
5138
5139  CoinAbcGatherFrom(abcCost_,costBasic_,abcPivotVariable_,numberRows_);
5140  return 0;
5141}
5142#if ABC_NORMAL_DEBUG>0
5143extern int cilk_conflict;
5144#endif
5145int 
5146AbcSimplexDual::noPivotRow()
5147{
5148#if ABC_NORMAL_DEBUG>3
5149  if (handler_->logLevel() & 32)
5150    printf("** no row pivot\n");
5151#endif
5152  int numberPivots = abcFactorization_->pivots();
5153  bool specialCase;
5154  int useNumberFake;
5155  int returnCode = 0;
5156  if (numberPivots <= CoinMax(dontFactorizePivots_, 20) &&
5157      (specialOptions_ & 2048) != 0 && currentDualBound_ >= 1.0e8) {
5158    specialCase = true;
5159    // as dual bound high - should be okay
5160    useNumberFake = 0;
5161  } else {
5162    specialCase = false;
5163    useNumberFake = numberAtFakeBound();
5164  }
5165  if (!numberPivots || specialCase) {
5166    // may have crept through - so may be optimal
5167    // check any flagged variables
5168    if (numberFlagged_ && numberPivots) {
5169      // try factorization
5170      returnCode = -2;
5171    }
5172   
5173    if (useNumberFake || numberDualInfeasibilities_) {
5174      // may be dual infeasible
5175      if ((specialOptions_ & 1024) == 0)
5176        problemStatus_ = -5;
5177      else if (!useNumberFake && numberPrimalInfeasibilities_
5178               && !numberPivots)
5179        problemStatus_ = 1;
5180    } else {
5181      if (numberFlagged_) {
5182#if ABC_NORMAL_DEBUG>3
5183        std::cout << "Flagged variables at end - infeasible?" << std::endl;
5184        printf("Probably infeasible - pivot was %g\n", alpha_);
5185#endif
5186#if ABC_NORMAL_DEBUG>3
5187        abort();
5188#endif
5189        problemStatus_ = -5;
5190      } else {
5191        problemStatus_ = 0;
5192#ifndef CLP_CHECK_NUMBER_PIVOTS
5193#define CLP_CHECK_NUMBER_PIVOTS 10
5194#endif
5195#if CLP_CHECK_NUMBER_PIVOTS < 20
5196        if (numberPivots > CLP_CHECK_NUMBER_PIVOTS) {
5197          gutsOfSolution(2);
5198          rawObjectiveValue_ +=sumNonBasicCosts_;
5199          setClpSimplexObjectiveValue();
5200          if (numberPrimalInfeasibilities_) {
5201            problemStatus_ = -1;
5202            returnCode = -2;
5203          }
5204        }
5205#endif
5206        if (!problemStatus_) {
5207          // make it look OK
5208          numberPrimalInfeasibilities_ = 0;
5209          sumPrimalInfeasibilities_ = 0.0;
5210          numberDualInfeasibilities_ = 0;
5211          sumDualInfeasibilities_ = 0.0;
5212          if (numberChanged_) {
5213            numberChanged_ = 0; // Number of variables with changed costs
5214            CoinAbcMemcpy(abcCost_,abcCost_+maximumNumberTotal_,numberTotal_);
5215            CoinAbcGatherFrom(abcCost_,costBasic_,abcPivotVariable_,numberRows_);
5216                  // make sure duals are current
5217            gutsOfSolution(1);
5218            setClpSimplexObjectiveValue();
5219            abcProgress_.modifyObjective(-COIN_DBL_MAX);
5220            if (numberDualInfeasibilities_) {
5221              problemStatus_ = 10; // was -3;
5222              numberTimesOptimal_++;
5223            } else {
5224              computeObjectiveValue(true);
5225            }
5226          } else if (numberPivots) {
5227            computeObjectiveValue(true);
5228          }
5229          sumPrimalInfeasibilities_ = 0.0;
5230        }
5231        if ((specialOptions_&(1024 + 16384)) != 0 && !problemStatus_&&abcFactorization_->pivots()) {
5232          CoinIndexedVector * arrayVector = &usefulArray_[arrayForFtran_];
5233          double * rhs = arrayVector->denseVector();
5234          arrayVector->clear();
5235          CoinAbcScatterTo(solutionBasic_,abcSolution_,abcPivotVariable_,numberRows_);
5236          abcMatrix_->timesModifyExcludingSlacks(-1.0, abcSolution_,rhs);
5237          //#define CHECK_ACCURACY
5238#ifdef CHECK_ACCURACY
5239          bool bad = false;
5240#endif
5241          bool bad2 = false;
5242          int i;
5243          for ( i = 0; i < numberRows_; i++) {
5244            if (rhs[i] < abcLower_[i] - primalTolerance_ ||
5245                rhs[i] > abcUpper_[i] + primalTolerance_) {
5246              bad2 = true;
5247#ifdef CHECK_ACCURACY
5248              printf("row %d out of bounds %g, %g correct %g bad %g\n", i,
5249                     abcLower_[i], abcUpper_[i],
5250                     rhs[i], abcSolution_[i]);
5251#endif
5252            } else if (fabs(rhs[i] - abcSolution_[i]) > 1.0e-3) {
5253#ifdef CHECK_ACCURACY
5254              bad = true;
5255              printf("row %d correct %g bad %g\n", i, rhs[i], abcSolution_[i]);
5256#endif
5257            }
5258            rhs[i] = 0.0;
5259          }
5260          for (int i = 0; i < numberTotal_; i++) {
5261            if (abcSolution_[i] < abcLower_[i] - primalTolerance_ ||
5262                abcSolution_[i] > abcUpper_[i] + primalTolerance_) {
5263              bad2 = true;
5264#ifdef CHECK_ACCURACY
5265              printf("column %d out of bounds %g, %g value %g\n", i,
5266                     abcLower_[i], abcUpper_[i],
5267                     abcSolution_[i]);
5268#endif
5269            }
5270          }
5271          if (bad2) {
5272            problemStatus_ = -3;
5273            returnCode = -2;
5274            // Force to re-factorize early next time
5275            int numberPivots = abcFactorization_->pivots();
5276            if (forceFactorization_<0)
5277              forceFactorization_=100000;
5278            forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
5279          }
5280        }
5281      }
5282    }
5283  } else {
5284    problemStatus_ = -3;
5285    returnCode = -2;
5286    // Force to re-factorize early next time
5287    int numberPivots = abcFactorization_->pivots();
5288    if (forceFactorization_<0)
5289      forceFactorization_=100000;
5290    forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
5291  }
5292  return returnCode;
5293}
5294void 
5295AbcSimplexDual::dualPivotColumn()
5296{
5297  double saveAcceptable=acceptablePivot_;
5298  if (largestPrimalError_>1.0e-5||largestDualError_>1.0e-5) {
5299    //if ((largestPrimalError_>1.0e-5||largestDualError_>1.0e-5)&&false) {
5300    if (!abcFactorization_->pivots())
5301      acceptablePivot_*=1.0e2;
5302    else if (abcFactorization_->pivots()<5)
5303      acceptablePivot_*=1.0e1;
5304  }
5305  dualColumn2();
5306  if (sequenceIn_<0&&acceptablePivot_>saveAcceptable) {
5307    acceptablePivot_=saveAcceptable;
5308#if ABC_NORMAL_DEBUG>0
5309    printf("No pivot!!!\n");
5310#endif
5311    // try again
5312    dualColumn1(true);
5313    dualColumn2();
5314  } else {
5315    acceptablePivot_=saveAcceptable;
5316  }
5317  if ((stateOfProblem_&VALUES_PASS)!=0) {
5318    // check no flips
5319    assert (!usefulArray_[arrayForFlipBounds_].getNumElements());
5320    // If no pivot - then try other way if plausible
5321    if (sequenceIn_<0) {
5322      // could use fake primal basic values
5323      if ((directionOut_<0&&lowerOut_>-1.0e30)||
5324          (directionOut_>0&&upperOut_<1.30)) {
5325        directionOut_=-directionOut_;
5326        if (directionOut_<0&&abcDj_[sequenceOut_]<0.0)
5327          upperTheta_=-abcDj_[sequenceOut_];
5328        else if (directionOut_>0&&abcDj_[sequenceOut_]>0.0)
5329          upperTheta_=abcDj_[sequenceOut_];
5330        else
5331          upperTheta_=1.0e30;
5332        CoinPartitionedVector & tableauRow = usefulArray_[arrayForTableauRow_];
5333        CoinPartitionedVector & candidateList = usefulArray_[arrayForDualColumn_];
5334        assert (!candidateList.getNumElements());
5335        tableauRow.compact();
5336        candidateList.compact();
5337        // redo candidate list
5338        int *  COIN_RESTRICT index = tableauRow.getIndices();
5339        double *  COIN_RESTRICT array = tableauRow.denseVector();
5340        double *  COIN_RESTRICT arrayCandidate=candidateList.denseVector();
5341        int *  COIN_RESTRICT indexCandidate = candidateList.getIndices();
5342        const double *  COIN_RESTRICT abcDj = abcDj_;
5343        const unsigned char *  COIN_RESTRICT internalStatus = internalStatus_;
5344        // do first pass to get possibles
5345        double bestPossible = 0.0;
5346        // We can also see if infeasible or pivoting on free
5347        double tentativeTheta = 1.0e25; // try with this much smaller as guess
5348        double acceptablePivot = currentAcceptablePivot_;
5349        double dualT=-currentDualTolerance_;
5350        // fixed will have been taken out by now
5351        const double multiplier[] = { 1.0, -1.0};
5352        int freeSequence=-1;
5353        double freeAlpha=0.0;
5354        int numberNonZero=tableauRow.getNumElements();
5355        int numberRemaining=0;
5356        if (ordinaryVariables()) {
5357          for (int i = 0; i < numberNonZero; i++) {
5358            int iSequence = index[i];
5359            double tableauValue=array[i];
5360            unsigned char iStatus=internalStatus[iSequence]&7;
5361            double mult = multiplier[iStatus];
5362            double alpha = tableauValue * mult;
5363            double oldValue = abcDj[iSequence] * mult;
5364            double value = oldValue - tentativeTheta * alpha;
5365            if (value < dualT) {
5366              bestPossible = CoinMax(bestPossible, alpha);
5367              value = oldValue - upperTheta_ * alpha;
5368              if (value < dualT && alpha >= acceptablePivot) {
5369                upperTheta_ = (oldValue - dualT) / alpha;
5370              }
5371              // add to list
5372              arrayCandidate[numberRemaining] = alpha;
5373              indexCandidate[numberRemaining++] = iSequence;
5374            }
5375          }
5376        } else {
5377          double badFree = 0.0;
5378          freeAlpha = currentAcceptablePivot_;
5379          double currentDualTolerance = currentDualTolerance_;
5380          for (int i = 0; i < numberNonZero; i++) {
5381            int iSequence = index[i];
5382            double tableauValue=array[i];
5383            unsigned char iStatus=internalStatus[iSequence]&7;
5384            if ((iStatus&4)==0) {
5385              double mult = multiplier[iStatus];
5386              double alpha = tableauValue * mult;
5387              double oldValue = abcDj[iSequence] * mult;
5388              double value = oldValue - tentativeTheta * alpha;
5389              if (value < dualT) {
5390                bestPossible = CoinMax(bestPossible, alpha);
5391                value = oldValue - upperTheta_ * alpha;
5392                if (value < dualT && alpha >= acceptablePivot) {
5393                  upperTheta_ = (oldValue - dualT) / alpha;
5394                }
5395                // add to list
5396                arrayCandidate[numberRemaining] = alpha;
5397                indexCandidate[numberRemaining++] = iSequence;
5398              }
5399            } else {
5400              bool keep;
5401              bestPossible = CoinMax(bestPossible, fabs(tableauValue));
5402              double oldValue = abcDj[iSequence];
5403              // If free has to be very large - should come in via dualRow
5404              //if (getInternalStatus(iSequence+addSequence)==isFree&&fabs(tableauValue)<1.0e-3)
5405              //break;
5406              if (oldValue > currentDualTolerance) {
5407                keep = true;
5408              } else if (oldValue < -currentDualTolerance) {
5409                keep = true;
5410              } else {
5411                if (fabs(tableauValue) > CoinMax(10.0 * acceptablePivot, 1.0e-5)) {
5412                  keep = true;
5413                } else {
5414                  keep = false;
5415                  badFree = CoinMax(badFree, fabs(tableauValue));
5416                }
5417              }
5418              if (keep) {
5419                // free - choose largest
5420                if (fabs(tableauValue) > fabs(freeAlpha)) {
5421                  freeAlpha = tableauValue;
5422                  freeSequence = iSequence;
5423                }
5424              }
5425            }
5426          }
5427        }
5428        candidateList.setNumElements(numberRemaining);
5429        if (freeSequence>=0) {
5430          sequenceIn_=freeSequence;
5431        } else if (fabs(upperTheta_-fabs(abcDj_[sequenceOut_]))<dualTolerance_) {
5432          // can just move
5433          stateOfIteration_=-1;
5434          return;
5435        }
5436        dualColumn2();
5437      }
5438    }
5439    // redo dualOut_
5440    if (directionOut_<0) {
5441      dualOut_ = valueOut_ - upperOut_;
5442    } else {
5443      dualOut_ = lowerOut_ - valueOut_;
5444    }
5445  }
5446  usefulArray_[arrayForDualColumn_].clear();
5447  //clearArrays(arrayForDualColumn_);
5448  if (sequenceIn_ >= 0&&upperTheta_==COIN_DBL_MAX) {
5449    lowerIn_ = abcLower_[sequenceIn_];
5450    upperIn_ = abcUpper_[sequenceIn_];
5451    valueIn_ = abcSolution_[sequenceIn_];
5452    dualIn_ = abcDj_[sequenceIn_];
5453    if (valueIn_<lowerIn_+primalTolerance_) {
5454      directionIn_ = 1;
5455    } else if (valueIn_>upperIn_-primalTolerance_) {
5456      directionIn_ = -1;
5457    } else {
5458      if (alpha_ < 0.0) {
5459        // as if from upper bound
5460        directionIn_ = -1;
5461        //assert(upperIn_ == valueIn_);
5462      } else {
5463        // as if from lower bound
5464        directionIn_ = 1;
5465        //assert(lowerIn_ == valueIn_);
5466      }
5467    }
5468  }
5469#if ABC_DEBUG
5470  checkArrays(4);
5471#endif
5472#if CAN_HAVE_ZERO_OBJ>1
5473  if ((specialOptions_&2097152)!=0) 
5474    theta_=0.0;
5475#endif
5476  movement_=0.0;
5477  objectiveChange_=0.0;
5478  /*
5479    0 - take iteration
5480    1 - don't take but continue
5481    2 - don't take and break
5482  */
5483  btranAlpha_ = -alpha_ * directionOut_; // for check
5484}
5485void
5486AbcSimplexDual::getTableauColumnPart2()
5487{
5488#if ABC_PARALLEL
5489  if (parallelMode_!=0) {
5490    abcFactorization_->updateColumnFTPart2(usefulArray_[arrayForFtran_]);
5491    // pivot element
5492    alpha_ = usefulArray_[arrayForFtran_].denseVector()[pivotRow_];
5493  }
5494#endif
5495}
5496void 
5497AbcSimplexDual::replaceColumnPart3()
5498{
5499  abcFactorization_->replaceColumnPart3(this,
5500                                        &usefulArray_[arrayForReplaceColumn_],
5501                                        &usefulArray_[arrayForFtran_],
5502                                        lastPivotRow_,
5503                                        ftAlpha_?ftAlpha_:alpha_);
5504}
5505void
5506AbcSimplexDual::checkReplacePart1()
5507{
5508  //abcFactorization_->checkReplacePart1a(&usefulArray_[arrayForReplaceColumn_],pivotRow_);
5509  //usefulArray_[arrayForReplaceColumn_].print();
5510  ftAlpha_=abcFactorization_->checkReplacePart1(&usefulArray_[arrayForReplaceColumn_],pivotRow_);
5511}
5512#if 0
5513void
5514AbcSimplexDual::checkReplacePart1a()
5515{
5516  abcFactorization_->checkReplacePart1a(&usefulArray_[arrayForReplaceColumn_],pivotRow_);
5517}
5518void
5519AbcSimplexDual::checkReplacePart1b()
5520{
5521  //usefulArray_[arrayForReplaceColumn_].print();
5522  ftAlpha_=abcFactorization_->checkReplacePart1b(&usefulArray_[arrayForReplaceColumn_],pivotRow_);
5523}
5524#endif
5525/*
5526  returns
5527  0 - OK
5528  1 - take iteration then re-factorize
5529  2 - flag something and skip
5530  3 - break and re-factorize
5531  5 - take iteration then re-factorize because of memory
5532 */
5533int 
5534AbcSimplexDual::checkReplace()
5535{
5536  int returnCode=0;
5537  if (!stateOfIteration_) {
5538    int saveStatus=problemStatus_;
5539    // check update
5540#if MOVE_REPLACE_PART1A < 0
5541    checkReplacePart1();
5542#endif
5543    int updateStatus=abcFactorization_->checkReplacePart2(pivotRow_,btranAlpha_,alpha_,ftAlpha_);
5544#if ABC_DEBUG
5545    checkArrays();
5546#endif
5547    // If looks like bad pivot - refactorize
5548    if (fabs(dualOut_) > 1.0e50)
5549      updateStatus = 2;
5550    // if no pivots, bad update but reasonable alpha - take and invert
5551    if (updateStatus == 2 &&
5552        !abcFactorization_->pivots() && fabs(alpha_) > 1.0e-5)
5553      updateStatus = 4;
5554    if (updateStatus == 1 || updateStatus == 4) {
5555      // slight error
5556      if (abcFactorization_->pivots() > 5 || updateStatus == 4) {
5557        problemStatus_ = -2; // factorize now
5558        returnCode = -3;
5559      }
5560    } else if (updateStatus == 2) {
5561      // major error
5562      // later we may need to unwind more e.g. fake bounds
5563      if (abcFactorization_->pivots() &&
5564          ((moreSpecialOptions_ & 16) == 0 || abcFactorization_->pivots() > 4)) {
5565        problemStatus_ = -2; // factorize now
5566        returnCode = -2;
5567        moreSpecialOptions_ |= 16;
5568        stateOfIteration_=2;
5569      } else {
5570        assert (sequenceOut_>=0);
5571        // need to reject something
5572        char x = isColumn(sequenceOut_) ? 'C' : 'R';
5573        handler_->message(CLP_SIMPLEX_FLAG, messages_)
5574          << x << sequenceWithin(sequenceOut_)
5575          << CoinMessageEol;
5576#if ABC_NORMAL_DEBUG>0
5577        printf("BAD flag b %g pivotRow %d (seqout %d) sequenceIn %d\n", alpha_,pivotRow_,
5578               sequenceOut_,sequenceIn_);
5579#endif
5580        setFlagged(sequenceOut_);
5581        // so can't happen immediately again
5582        sequenceOut_=-1;
5583        //abcProgress_.clearBadTimes();
5584        lastBadIteration_ = numberIterations_; // say be more cautious
5585        //clearArrays(-1);
5586        //if(primalSolutionUpdated)
5587        //gutsOfSolution(2);
5588        stateOfIteration_=1;
5589      }
5590    } else if (updateStatus == 3) {
5591      // out of memory
5592      // increase space if not many iterations
5593      if (abcFactorization_->pivots() <
5594          0.5 * abcFactorization_->maximumPivots() &&
5595          abcFactorization_->pivots() < 200)
5596        abcFactorization_->areaFactor(
5597                                      abcFactorization_->areaFactor() * 1.1);
5598      problemStatus_ = -2; // factorize now
5599    } else if (updateStatus == 5) {
5600      problemStatus_ = -2; // factorize now
5601<