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

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

minor changes to implement Aboca

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