# source:trunk/Clp/src/AbcSimplexDual.cpp@1879

Last change on this file since 1879 was 1879, checked in by forrest, 9 years ago

fix a few problems with Aboca

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