source: trunk/Clp/src/ClpSimplexDual.cpp @ 1394

Last change on this file since 1394 was 1394, checked in by forrest, 10 years ago

allow infeasibilityray more and prepare

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