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

Last change on this file since 1294 was 1294, checked in by forrest, 11 years ago

make numberPrimalInfeasibilities_ more plausible

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