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

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

modify asserts etc for difficult problem

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