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

Last change on this file since 1344 was 1344, checked in by forrest, 12 years ago

changes to simplex and lots of stuff and start Mumps cholesky

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 210.3 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              abort();
2722              numberFake_++;
2723            } else {
2724              setFakeBound(iSequence,ClpSimplexDual::lowerFake);
2725              numberFake_++;
2726            }
2727          } else {
2728            if (newUpperValue < upperValue) {
2729              setFakeBound(iSequence,ClpSimplexDual::upperFake);
2730              numberFake_++;
2731            }
2732          }
2733          if (status==atUpperBound)
2734            solution_[iSequence] = newUpperValue;
2735          else
2736            solution_[iSequence] = newLowerValue;
2737          double movement = solution_[iSequence] - value;
2738          if (movement&&outputArray) {
2739            if (iSequence>=numberColumns_) {
2740              outputArray->quickAdd(iSequence,-movement);
2741              changeCost += movement*cost_[iSequence];
2742            } else {
2743              matrix_->add(this,outputArray,iSequence,movement);
2744              changeCost += movement*cost_[iSequence];
2745            }
2746          }
2747        }
2748      }
2749      dualBound_ = newBound;
2750    } else {
2751      numberInfeasibilities=-1;
2752    }
2753    return numberInfeasibilities;
2754  } else if (initialize==1||initialize==3) {
2755    int iSequence;
2756    if (initialize==3) {
2757      for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
2758        setFakeBound(iSequence,ClpSimplexDual::noFake);
2759      }
2760    }
2761    double testBound = 0.999999*dualBound_;
2762    for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
2763      Status status = getStatus(iSequence);
2764      if (status==atUpperBound||
2765          status==atLowerBound) {
2766        double lowerValue=lower_[iSequence];
2767        double upperValue=upper_[iSequence];
2768        double value = solution_[iSequence];
2769        if (lowerValue>-largeValue_||upperValue<largeValue_) {
2770          if (true||lowerValue-value>-0.5*dualBound_||
2771              upperValue-value<0.5*dualBound_) {
2772            if (fabs(lowerValue-value)<=fabs(upperValue-value)) {
2773              if (upperValue > lowerValue + testBound) {
2774                if (getFakeBound(iSequence)==ClpSimplexDual::noFake)
2775                  numberFake_++;
2776                upper_[iSequence]=lowerValue+dualBound_;
2777                setFakeBound(iSequence,ClpSimplexDual::upperFake);
2778              }
2779            } else {
2780              if (lowerValue < upperValue - testBound) {
2781                if (getFakeBound(iSequence)==ClpSimplexDual::noFake)
2782                  numberFake_++;
2783                lower_[iSequence]=upperValue-dualBound_;
2784                setFakeBound(iSequence,ClpSimplexDual::lowerFake);
2785              }
2786            }
2787          } else {
2788            if (getFakeBound(iSequence)==ClpSimplexDual::noFake)
2789              numberFake_++;
2790            lower_[iSequence]=-0.5*dualBound_;
2791            upper_[iSequence]= 0.5*dualBound_;
2792            setFakeBound(iSequence,ClpSimplexDual::bothFake);
2793            abort();
2794          }
2795          if (status==atUpperBound)
2796            solution_[iSequence]=upper_[iSequence];
2797          else
2798            solution_[iSequence]=lower_[iSequence];
2799        } else {
2800          // set non basic free variables to fake bounds
2801          // I don't think we should ever get here
2802          CoinAssert(!("should not be here"));
2803          lower_[iSequence]=-0.5*dualBound_;
2804          upper_[iSequence]= 0.5*dualBound_;
2805          setFakeBound(iSequence,ClpSimplexDual::bothFake);
2806          abort();
2807          numberFake_++;
2808          setStatus(iSequence,atUpperBound);
2809          solution_[iSequence]=0.5*dualBound_;
2810        }
2811      } else if (status==basic) {
2812        // make sure not at fake bound and bounds correct
2813        setFakeBound(iSequence,ClpSimplexDual::noFake);
2814        double gap = upper_[iSequence]-lower_[iSequence];
2815        if (gap>0.5*dualBound_&&gap<2.0*dualBound_) {
2816          if (iSequence<numberColumns_) {
2817            if (columnScale_) {
2818              double multiplier = rhsScale_*inverseColumnScale_[iSequence];
2819              // lower
2820              double value = columnLower_[iSequence];
2821              if (value>-1.0e30) {
2822                value *= multiplier;
2823              }
2824              lower_[iSequence]=value;
2825              // upper
2826              value = columnUpper_[iSequence];
2827              if (value<1.0e30) {
2828                value *= multiplier;
2829              }
2830              upper_[iSequence]=value;
2831            } else {
2832              lower_[iSequence]=columnLower_[iSequence];;
2833              upper_[iSequence]=columnUpper_[iSequence];;
2834            }
2835          } else {
2836            int iRow = iSequence-numberColumns_;
2837            if (rowScale_) {
2838              // lower
2839              double multiplier = rhsScale_*rowScale_[iRow];
2840              double value = rowLower_[iRow];
2841              if (value>-1.0e30) {
2842                value *= multiplier;
2843              }
2844              lower_[iSequence]=value;
2845              // upper
2846              value = rowUpper_[iRow];
2847              if (value<1.0e30) {
2848                value *= multiplier;
2849              }
2850              upper_[iSequence]=value;
2851            } else {
2852              lower_[iSequence]=rowLower_[iRow];;
2853              upper_[iSequence]=rowUpper_[iRow];;
2854            }
2855          }
2856        }
2857      }
2858    }
2859
2860    return 1;
2861  } else {
2862    // just reset changed ones
2863    if (columnScale_) {
2864      int iSequence;
2865      for (iSequence=0;iSequence<numberColumns_;iSequence++) {
2866        FakeBound fakeStatus = getFakeBound(iSequence);
2867        if (fakeStatus!=noFake) {
2868          if ((static_cast<int> (fakeStatus)&1)!=0) {
2869            // lower
2870            double value = columnLower_[iSequence];
2871            if (value>-1.0e30) {
2872              double multiplier = rhsScale_*inverseColumnScale_[iSequence];
2873              value *= multiplier;
2874            }
2875            columnLowerWork_[iSequence]=value;
2876          }
2877          if ((static_cast<int> (fakeStatus)&2)!=0) {
2878            // upper
2879            double value = columnUpper_[iSequence];
2880            if (value<1.0e30) {
2881              double multiplier = rhsScale_*inverseColumnScale_[iSequence];
2882              value *= multiplier;
2883            }
2884            columnUpperWork_[iSequence]=value;
2885          }
2886        }
2887      }
2888      for (iSequence=0;iSequence<numberRows_;iSequence++) {
2889        FakeBound fakeStatus = getFakeBound(iSequence+numberColumns_);
2890        if (fakeStatus!=noFake) {
2891          if ((static_cast<int> (fakeStatus)&1)!=0) {
2892            // lower
2893            double value = rowLower_[iSequence];
2894            if (value>-1.0e30) {
2895              double multiplier = rhsScale_*rowScale_[iSequence];
2896              value *= multiplier;
2897            }
2898            rowLowerWork_[iSequence]=value;
2899          }
2900          if ((static_cast<int> (fakeStatus)&2)!=0) {
2901            // upper
2902            double value = rowUpper_[iSequence];
2903            if (value<1.0e30) {
2904              double multiplier = rhsScale_*rowScale_[iSequence];
2905              value *= multiplier;
2906            }
2907            rowUpperWork_[iSequence]=value;
2908          }
2909        }
2910      }
2911    } else {
2912      int iSequence;
2913      for (iSequence=0;iSequence<numberColumns_;iSequence++) {
2914        FakeBound fakeStatus = getFakeBound(iSequence);
2915        if ((static_cast<int> (fakeStatus)&1)!=0) {
2916          // lower
2917          columnLowerWork_[iSequence]=columnLower_[iSequence];
2918        }
2919        if ((static_cast<int> (fakeStatus)&2)!=0) {
2920          // upper
2921          columnUpperWork_[iSequence]=columnUpper_[iSequence];
2922        }
2923      }
2924      for (iSequence=0;iSequence<numberRows_;iSequence++) {
2925        FakeBound fakeStatus = getFakeBound(iSequence+numberColumns_);
2926        if ((static_cast<int> (fakeStatus)&1)!=0) {
2927          // lower
2928          rowLowerWork_[iSequence]=rowLower_[iSequence];
2929        }
2930        if ((static_cast<int> (fakeStatus)&2)!=0) {
2931          // upper
2932          rowUpperWork_[iSequence]=rowUpper_[iSequence];
2933        }
2934      }
2935    }
2936    return 0;
2937  }
2938}
2939int
2940ClpSimplexDual::dualColumn0(const CoinIndexedVector * rowArray,
2941                           const CoinIndexedVector * columnArray,
2942                           CoinIndexedVector * spareArray,
2943                           double acceptablePivot,
2944                           double & upperReturn, double &bestReturn,double & badFree)
2945{
2946  // do first pass to get possibles
2947  double * spare = spareArray->denseVector();
2948  int * index = spareArray->getIndices();
2949  const double * work;
2950  int number;
2951  const int * which;
2952  const double * reducedCost;
2953  // We can also see if infeasible or pivoting on free
2954  double tentativeTheta = 1.0e25;
2955  double upperTheta = 1.0e31;
2956  double freePivot = acceptablePivot;
2957  double bestPossible=0.0;
2958  int numberRemaining=0;
2959  int i;
2960  badFree=0.0;
2961  for (int iSection=0;iSection<2;iSection++) {
2962
2963    int addSequence;
2964
2965    if (!iSection) {
2966      work = rowArray->denseVector();
2967      number = rowArray->getNumElements();
2968      which = rowArray->getIndices();
2969      reducedCost = rowReducedCost_;
2970      addSequence = numberColumns_;
2971    } else {
2972      work = columnArray->denseVector();
2973      number = columnArray->getNumElements();
2974      which = columnArray->getIndices();
2975      reducedCost = reducedCostWork_;
2976      addSequence = 0;
2977    }
2978   
2979    for (i=0;i<number;i++) {
2980      int iSequence = which[i];
2981      double alpha;
2982      double oldValue;
2983      double value;
2984      bool keep;
2985
2986      switch(getStatus(iSequence+addSequence)) {
2987         
2988      case basic:
2989      case ClpSimplex::isFixed:
2990        break;
2991      case isFree:
2992      case superBasic:
2993        alpha = work[i];
2994        bestPossible = CoinMax(bestPossible,fabs(alpha));
2995        oldValue = reducedCost[iSequence];
2996        // If free has to be very large - should come in via dualRow
2997        //if (getStatus(iSequence+addSequence)==isFree&&fabs(alpha)<1.0e-3)
2998        //break;
2999        if (oldValue>dualTolerance_) {
3000          keep = true;
3001        } else if (oldValue<-dualTolerance_) {
3002          keep = true;
3003        } else {
3004          if (fabs(alpha)>CoinMax(10.0*acceptablePivot,1.0e-5)) {
3005            keep = true;
3006          } else {
3007            keep=false;
3008            badFree=CoinMax(badFree,fabs(alpha));
3009          }
3010        }
3011        if (keep) {
3012          // free - choose largest
3013          if (fabs(alpha)>freePivot) {
3014            freePivot=fabs(alpha);
3015            sequenceIn_ = iSequence + addSequence;
3016            theta_=oldValue/alpha;
3017            alpha_=alpha;
3018          }
3019        }
3020        break;
3021      case atUpperBound:
3022        alpha = work[i];
3023        oldValue = reducedCost[iSequence];
3024        value = oldValue-tentativeTheta*alpha;
3025        //assert (oldValue<=dualTolerance_*1.0001);
3026        if (value>dualTolerance_) {
3027          bestPossible = CoinMax(bestPossible,-alpha);
3028          value = oldValue-upperTheta*alpha;
3029          if (value>dualTolerance_ && -alpha>=acceptablePivot) 
3030            upperTheta = (oldValue-dualTolerance_)/alpha;
3031          // add to list
3032          spare[numberRemaining]=alpha;
3033          index[numberRemaining++]=iSequence+addSequence;
3034        }
3035        break;
3036      case atLowerBound:
3037        alpha = work[i];
3038        oldValue = reducedCost[iSequence];
3039        value = oldValue-tentativeTheta*alpha;
3040        //assert (oldValue>=-dualTolerance_*1.0001);
3041        if (value<-dualTolerance_) {
3042          bestPossible = CoinMax(bestPossible,alpha);
3043          value = oldValue-upperTheta*alpha;
3044          if (value<-dualTolerance_ && alpha>=acceptablePivot) 
3045            upperTheta = (oldValue+dualTolerance_)/alpha;
3046          // add to list
3047          spare[numberRemaining]=alpha;
3048          index[numberRemaining++]=iSequence+addSequence;
3049        }
3050        break;
3051      }
3052    }
3053  }
3054  upperReturn = upperTheta;
3055  bestReturn = bestPossible;
3056  return numberRemaining;
3057}
3058/*
3059   Row array has row part of pivot row (as duals so sign may be switched)
3060   Column array has column part.
3061   This chooses pivot column.
3062   Spare array will be needed when we start getting clever.
3063   We will check for basic so spare array will never overflow.
3064   If necessary will modify costs
3065*/
3066double
3067ClpSimplexDual::dualColumn(CoinIndexedVector * rowArray,
3068                           CoinIndexedVector * columnArray,
3069                           CoinIndexedVector * spareArray,
3070                           CoinIndexedVector * spareArray2,
3071                           double acceptablePivot,
3072                           CoinBigIndex * dubiousWeights)
3073{
3074  int numberPossiblySwapped=0;
3075  int numberRemaining=0;
3076 
3077  double totalThru=0.0; // for when variables flip
3078  //double saveAcceptable=acceptablePivot;
3079  //acceptablePivot=1.0e-9;
3080
3081  double bestEverPivot=acceptablePivot;
3082  int lastSequence = -1;
3083  double lastPivot=0.0;
3084  double upperTheta;
3085  double newTolerance = dualTolerance_;
3086  //newTolerance = dualTolerance_+1.0e-6*dblParam_[ClpDualTolerance];
3087  // will we need to increase tolerance
3088  bool thisIncrease=false;
3089  // If we think we need to modify costs (not if something from broad sweep)
3090  bool modifyCosts=false;
3091  // Increase in objective due to swapping bounds (may be negative)
3092  double increaseInObjective=0.0;
3093
3094  // use spareArrays to put ones looked at in
3095  // we are going to flip flop between
3096  int iFlip = 0;
3097  // Possible list of pivots
3098  int interesting[2];
3099  // where possible swapped ones are
3100  int swapped[2];
3101  // for zeroing out arrays after
3102  int marker[2][2];
3103  // pivot elements
3104  double * array[2], * spare, * spare2;
3105  // indices
3106  int * indices[2], * index, * index2;
3107  spareArray2->clear();
3108  array[0] = spareArray->denseVector();
3109  indices[0] = spareArray->getIndices();
3110  spare = array[0];
3111  index = indices[0];
3112  array[1] = spareArray2->denseVector();
3113  indices[1] = spareArray2->getIndices();
3114  int i;
3115
3116  // initialize lists
3117  for (i=0;i<2;i++) {
3118    interesting[i]=0;
3119    swapped[i]=numberColumns_;
3120    marker[i][0]=0;
3121    marker[i][1]=numberColumns_;
3122  }
3123  /*
3124    First we get a list of possible pivots.  We can also see if the
3125    problem looks infeasible or whether we want to pivot in free variable.
3126    This may make objective go backwards but can only happen a finite
3127    number of times and I do want free variables basic.
3128
3129    Then we flip back and forth.  At the start of each iteration
3130    interesting[iFlip] should have possible candidates and swapped[iFlip]
3131    will have pivots if we decide to take a previous pivot.
3132    At end of each iteration interesting[1-iFlip] should have
3133    candidates if we go through this theta and swapped[1-iFlip]
3134    pivots if we don't go through.
3135
3136    At first we increase theta and see what happens.  We start
3137    theta at a reasonable guess.  If in right area then we do bit by bit.
3138
3139   */
3140
3141  // do first pass to get possibles
3142  upperTheta = 1.0e31;
3143  double bestPossible=0.0;
3144  double badFree=0.0;
3145  alpha_=0.0;
3146  if (spareIntArray_[0]!=-1) {
3147    numberRemaining = dualColumn0(rowArray,columnArray,spareArray,
3148                                  acceptablePivot,upperTheta,bestPossible,badFree);
3149  } else {
3150    // already done
3151    numberRemaining = spareArray->getNumElements();
3152    spareArray->setNumElements(0);
3153    upperTheta = spareDoubleArray_[0];
3154    bestPossible = spareDoubleArray_[1];
3155    theta_ = spareDoubleArray_[2];
3156    alpha_ = spareDoubleArray_[3];
3157    sequenceIn_ = spareIntArray_[1];
3158  }
3159  // switch off
3160  spareIntArray_[0]=0;
3161  // We can also see if infeasible or pivoting on free
3162  double tentativeTheta = 1.0e25;
3163  interesting[0]=numberRemaining;
3164  marker[0][0] = numberRemaining;
3165
3166  if (!numberRemaining&&sequenceIn_<0)
3167    return 0.0; // Looks infeasible
3168
3169  // If sum of bad small pivots too much
3170#define MORE_CAREFUL
3171#ifdef MORE_CAREFUL
3172  bool badSumPivots=false;
3173#endif
3174  if (sequenceIn_>=0) {
3175    // free variable - always choose
3176  } else {
3177
3178    theta_=1.0e50;
3179    // now flip flop between spare arrays until reasonable theta
3180    tentativeTheta = CoinMax(10.0*upperTheta,1.0e-7);
3181   
3182    // loops increasing tentative theta until can't go through
3183   
3184    while (tentativeTheta < 1.0e22) {
3185      double thruThis = 0.0;
3186     
3187      double bestPivot=acceptablePivot;
3188      int bestSequence=-1;
3189     
3190      numberPossiblySwapped = numberColumns_;
3191      numberRemaining = 0;
3192
3193      upperTheta = 1.0e50;
3194
3195      spare = array[iFlip];
3196      index = indices[iFlip];
3197      spare2 = array[1-iFlip];
3198      index2 = indices[1-iFlip];
3199
3200      // try 3 different ways
3201      // 1 bias increase by ones with slightly wrong djs
3202      // 2 bias by all
3203      // 3 bias by all - tolerance
3204#define TRYBIAS 3
3205
3206
3207      double increaseInThis=0.0; //objective increase in this loop
3208     
3209      for (i=0;i<interesting[iFlip];i++) {
3210        int iSequence = index[i];
3211        double alpha = spare[i];
3212        double oldValue = dj_[iSequence];
3213        double value = oldValue-tentativeTheta*alpha;
3214
3215        if (alpha < 0.0) {
3216          //at upper bound
3217          if (value>newTolerance) {
3218            double range = upper_[iSequence] - lower_[iSequence];
3219            thruThis -= range*alpha;
3220#if TRYBIAS==1
3221            if (oldValue>0.0)
3222              increaseInThis -= oldValue*range;
3223#elif TRYBIAS==2
3224            increaseInThis -= oldValue*range;
3225#else
3226            increaseInThis -= (oldValue+dualTolerance_)*range;
3227#endif
3228            // goes on swapped list (also means candidates if too many)
3229            spare2[--numberPossiblySwapped]=alpha;
3230            index2[numberPossiblySwapped]=iSequence;
3231            if (fabs(alpha)>bestPivot) {
3232              bestPivot=fabs(alpha);
3233              bestSequence=numberPossiblySwapped;
3234            }
3235          } else {
3236            value = oldValue-upperTheta*alpha;
3237            if (value>newTolerance && -alpha>=acceptablePivot) 
3238              upperTheta = (oldValue-newTolerance)/alpha;
3239            spare2[numberRemaining]=alpha;
3240            index2[numberRemaining++]=iSequence;
3241          }
3242        } else {
3243          // at lower bound
3244          if (value<-newTolerance) {
3245            double range = upper_[iSequence] - lower_[iSequence];
3246            thruThis += range*alpha;
3247            //?? is this correct - and should we look at good ones
3248#if TRYBIAS==1
3249            if (oldValue<0.0)
3250              increaseInThis += oldValue*range;
3251#elif TRYBIAS==2
3252            increaseInThis += oldValue*range;
3253#else
3254            increaseInThis += (oldValue-dualTolerance_)*range;
3255#endif
3256            // goes on swapped list (also means candidates if too many)
3257            spare2[--numberPossiblySwapped]=alpha;
3258            index2[numberPossiblySwapped]=iSequence;
3259            if (fabs(alpha)>bestPivot) {
3260              bestPivot=fabs(alpha);
3261              bestSequence=numberPossiblySwapped;
3262            }
3263          } else {
3264            value = oldValue-upperTheta*alpha;
3265            if (value<-newTolerance && alpha>=acceptablePivot) 
3266              upperTheta = (oldValue+newTolerance)/alpha;
3267            spare2[numberRemaining]=alpha;
3268            index2[numberRemaining++]=iSequence;
3269          }
3270        }
3271      }
3272      swapped[1-iFlip]=numberPossiblySwapped;
3273      interesting[1-iFlip]=numberRemaining;
3274      marker[1-iFlip][0]= CoinMax(marker[1-iFlip][0],numberRemaining);
3275      marker[1-iFlip][1]= CoinMin(marker[1-iFlip][1],numberPossiblySwapped);
3276     
3277      if (totalThru+thruThis>=fabs(dualOut_)||
3278          increaseInObjective+increaseInThis<0.0) {
3279        // We should be pivoting in this batch
3280        // so compress down to this lot
3281        numberRemaining=0;
3282        for (i=numberColumns_-1;i>=swapped[1-iFlip];i--) {
3283          spare[numberRemaining]=spare2[i];
3284          index[numberRemaining++]=index2[i];
3285        }
3286        interesting[iFlip]=numberRemaining;
3287        int iTry;
3288#define MAXTRY 100
3289        // first get ratio with tolerance
3290        for (iTry=0;iTry<MAXTRY;iTry++) {
3291         
3292          upperTheta=1.0e50;
3293          numberPossiblySwapped = numberColumns_;
3294          numberRemaining = 0;
3295
3296          increaseInThis=0.0; //objective increase in this loop
3297
3298          thruThis=0.0;
3299
3300          spare = array[iFlip];
3301          index = indices[iFlip];
3302          spare2 = array[1-iFlip];
3303          index2 = indices[1-iFlip];
3304          for (i=0;i<interesting[iFlip];i++) {
3305            int iSequence=index[i];
3306            double alpha=spare[i];
3307            double oldValue = dj_[iSequence];
3308            double value = oldValue-upperTheta*alpha;
3309           
3310            if (alpha < 0.0) {
3311              //at upper bound
3312              if (value>newTolerance) {
3313                if (-alpha>=acceptablePivot) {
3314                  upperTheta = (oldValue-newTolerance)/alpha;
3315                }
3316              }
3317            } else {
3318              // at lower bound
3319              if (value<-newTolerance) {
3320                if (alpha>=acceptablePivot) {
3321                  upperTheta = (oldValue+newTolerance)/alpha;
3322                }
3323              }
3324            }
3325          }
3326          bestPivot=acceptablePivot;
3327          sequenceIn_=-1;
3328#ifdef DUBIOUS_WEIGHTS
3329          double bestWeight=COIN_DBL_MAX;
3330#endif
3331          double largestPivot=acceptablePivot;
3332          // now choose largest and sum all ones which will go through
3333          //printf("XX it %d number %d\n",numberIterations_,interesting[iFlip]);
3334  // Sum of bad small pivots
3335#ifdef MORE_CAREFUL
3336          double sumBadPivots=0.0;
3337          badSumPivots=false;
3338#endif
3339          // Make sure upperTheta will work (-O2 and above gives problems)
3340          upperTheta *= 1.0000000001;
3341          for (i=0;i<interesting[iFlip];i++) {
3342            int iSequence=index[i];
3343            double alpha=spare[i];
3344            double value = dj_[iSequence]-upperTheta*alpha;
3345            double badDj=0.0;
3346           
3347            bool addToSwapped=false;
3348           
3349            if (alpha < 0.0) {
3350              //at upper bound
3351              if (value>=0.0) { 
3352                addToSwapped=true;
3353#if TRYBIAS==1
3354                badDj = -CoinMax(dj_[iSequence],0.0);
3355#elif TRYBIAS==2
3356                badDj = -dj_[iSequence];
3357#else
3358                badDj = -dj_[iSequence]-dualTolerance_;
3359#endif
3360              }
3361            } else {
3362              // at lower bound
3363              if (value<=0.0) {
3364                addToSwapped=true;
3365#if TRYBIAS==1
3366                badDj = CoinMin(dj_[iSequence],0.0);
3367#elif TRYBIAS==2
3368                badDj = dj_[iSequence];
3369#else
3370                badDj = dj_[iSequence]-dualTolerance_;
3371#endif
3372              }
3373            }
3374            if (!addToSwapped) {
3375              // add to list of remaining
3376              spare2[numberRemaining]=alpha;
3377              index2[numberRemaining++]=iSequence;
3378            } else {
3379              // add to list of swapped
3380              spare2[--numberPossiblySwapped]=alpha;
3381              index2[numberPossiblySwapped]=iSequence;
3382              // select if largest pivot
3383              bool take=false;
3384              double absAlpha = fabs(alpha);
3385#ifdef DUBIOUS_WEIGHTS
3386              // User could do anything to break ties here
3387              double weight;
3388              if (dubiousWeights)
3389                weight=dubiousWeights[iSequence];
3390              else
3391                weight=1.0;
3392              weight += randomNumberGenerator_.randomDouble()*1.0e-2;
3393              if (absAlpha>2.0*bestPivot) {
3394                take=true;
3395              } else if (absAlpha>largestPivot) {
3396                // could multiply absAlpha and weight
3397                if (weight*bestPivot<bestWeight*absAlpha)
3398                  take=true;
3399              }
3400#else
3401              if (absAlpha>bestPivot) 
3402                take=true;
3403#endif
3404#ifdef MORE_CAREFUL
3405              if (absAlpha<acceptablePivot&&upperTheta<1.0e20) {
3406                if (alpha < 0.0) {
3407                  //at upper bound
3408                  if (value>dualTolerance_) {
3409                    double gap=upper_[iSequence]-lower_[iSequence];
3410                    if (gap<1.0e20)
3411                      sumBadPivots += value*gap; 
3412                    else
3413                      sumBadPivots += 1.0e20;
3414                    //printf("bad %d alpha %g dj at upper %g\n",
3415                    //     iSequence,alpha,value);
3416                  }
3417                } else {
3418                  //at lower bound
3419                  if (value<-dualTolerance_) {
3420                    double gap=upper_[iSequence]-lower_[iSequence];
3421                    if (gap<1.0e20)
3422                      sumBadPivots -= value*gap; 
3423                    else
3424                      sumBadPivots += 1.0e20;
3425                    //printf("bad %d alpha %g dj at lower %g\n",
3426                    //     iSequence,alpha,value);
3427                  }
3428                }
3429              }
3430#endif
3431#ifdef FORCE_FOLLOW
3432              if (iSequence==force_in) {
3433                printf("taking %d - alpha %g best %g\n",force_in,absAlpha,largestPivot);
3434                take=true;
3435              }
3436#endif
3437              if (take) {
3438                sequenceIn_ = numberPossiblySwapped;
3439                bestPivot =  absAlpha;
3440                theta_ = dj_[iSequence]/alpha;
3441                largestPivot = CoinMax(largestPivot,0.5*bestPivot);
3442#ifdef DUBIOUS_WEIGHTS
3443                bestWeight = weight;
3444#endif
3445                //printf(" taken seq %d alpha %g weight %d\n",
3446                //   iSequence,absAlpha,dubiousWeights[iSequence]);
3447              } else {
3448                //printf(" not taken seq %d alpha %g weight %d\n",
3449                //   iSequence,absAlpha,dubiousWeights[iSequence]);
3450              }
3451              double range = upper_[iSequence] - lower_[iSequence];
3452              thruThis += range*fabs(alpha);
3453              increaseInThis += badDj*range;
3454            }
3455          }
3456          marker[1-iFlip][0]= CoinMax(marker[1-iFlip][0],numberRemaining);
3457          marker[1-iFlip][1]= CoinMin(marker[1-iFlip][1],numberPossiblySwapped);
3458#ifdef MORE_CAREFUL
3459          // If we have done pivots and things look bad set alpha_ 0.0 to force factorization
3460          if (sumBadPivots>1.0e4) {
3461            if (handler_->logLevel()>1)
3462            printf("maybe forcing re-factorization - sum %g  %d pivots\n",sumBadPivots,
3463                   factorization_->pivots());
3464            if(factorization_->pivots()>3) {
3465              badSumPivots=true;
3466              break;
3467            }
3468          }
3469#endif
3470          swapped[1-iFlip]=numberPossiblySwapped;
3471          interesting[1-iFlip]=numberRemaining;
3472          // If we stop now this will be increase in objective (I think)
3473          double increase = (fabs(dualOut_)-totalThru)*theta_;
3474          increase += increaseInObjective;
3475          if (theta_<0.0)
3476            thruThis += fabs(dualOut_); // force using this one
3477          if (increaseInObjective<0.0&&increase<0.0&&lastSequence>=0) {
3478            // back
3479            // We may need to be more careful - we could do by
3480            // switch so we always do fine grained?
3481            bestPivot=0.0;
3482          } else {
3483            // add in
3484            totalThru += thruThis;
3485            increaseInObjective += increaseInThis;
3486          }
3487          if (bestPivot<0.1*bestEverPivot&&
3488              bestEverPivot>1.0e-6&&
3489              (bestPivot<1.0e-3||totalThru*2.0>fabs(dualOut_))) {
3490            // back to previous one
3491            sequenceIn_=lastSequence;
3492            // swap regions
3493            iFlip = 1-iFlip;
3494            break;
3495          } else if (sequenceIn_==-1&&upperTheta>largeValue_) {
3496            if (lastPivot>acceptablePivot) {
3497              // back to previous one
3498              sequenceIn_=lastSequence;
3499              // swap regions
3500              iFlip = 1-iFlip;
3501            } else {
3502              // can only get here if all pivots too small
3503            }
3504            break;
3505          } else if (totalThru>=fabs(dualOut_)) {
3506            modifyCosts=true; // fine grain - we can modify costs
3507            break; // no point trying another loop
3508          } else {
3509            lastSequence=sequenceIn_;
3510            if (bestPivot>bestEverPivot)
3511              bestEverPivot=bestPivot;
3512            iFlip = 1 -iFlip;
3513            modifyCosts=true; // fine grain - we can modify costs
3514          }
3515        }
3516        if (iTry==MAXTRY)
3517          iFlip = 1-iFlip; // flip back
3518        break;
3519      } else {
3520        // skip this lot
3521        if (bestPivot>1.0e-3||bestPivot>bestEverPivot) {
3522          bestEverPivot=bestPivot;
3523          lastSequence=bestSequence;
3524        } else {
3525          // keep old swapped
3526          CoinMemcpyN(array[iFlip]+swapped[iFlip],
3527                 numberColumns_-swapped[iFlip],array[1-iFlip]+swapped[iFlip]);
3528          CoinMemcpyN(indices[iFlip]+swapped[iFlip],
3529                 numberColumns_-swapped[iFlip],indices[1-iFlip]+swapped[iFlip]);
3530          marker[1-iFlip][1] = CoinMin(marker[1-iFlip][1],swapped[iFlip]);
3531          swapped[1-iFlip]=swapped[iFlip];
3532        }
3533        increaseInObjective += increaseInThis;
3534        iFlip = 1 - iFlip; // swap regions
3535        tentativeTheta = 2.0*upperTheta;
3536        totalThru += thruThis;
3537      }
3538    }
3539   
3540    // can get here without sequenceIn_ set but with lastSequence
3541    if (sequenceIn_<0&&lastSequence>=0) {
3542      // back to previous one
3543      sequenceIn_=lastSequence;
3544      // swap regions
3545      iFlip = 1-iFlip;
3546    }
3547   
3548#define MINIMUMTHETA 1.0e-18
3549    // Movement should be minimum for anti-degeneracy - unless
3550    // fixed variable out
3551    double minimumTheta;
3552    if (upperOut_>lowerOut_)
3553      minimumTheta=MINIMUMTHETA;
3554    else
3555      minimumTheta=0.0;
3556    if (sequenceIn_>=0) {
3557      // at this stage sequenceIn_ is just pointer into index array
3558      // flip just so we can use iFlip
3559      iFlip = 1 -iFlip;
3560      spare = array[iFlip];
3561      index = indices[iFlip];
3562      double oldValue;
3563      alpha_ = spare[sequenceIn_];
3564      sequenceIn_ = indices[iFlip][sequenceIn_];
3565      oldValue = dj_[sequenceIn_];
3566      theta_ = CoinMax(oldValue/alpha_,0.0);
3567      if (theta_<minimumTheta&&fabs(alpha_)<1.0e5&&1) {
3568        // can't pivot to zero
3569#if 0
3570        if (oldValue-minimumTheta*alpha_>=-dualTolerance_) {
3571          theta_=minimumTheta;
3572        } else if (oldValue-minimumTheta*alpha_>=-newTolerance) {
3573          theta_=minimumTheta;
3574          thisIncrease=true;
3575        } else {
3576          theta_=CoinMax((oldValue+newTolerance)/alpha_,0.0);
3577          thisIncrease=true;
3578        }
3579#else
3580        theta_=minimumTheta;
3581#endif
3582      }
3583      // may need to adjust costs so all dual feasible AND pivoted is exactly 0
3584      //int costOffset = numberRows_+numberColumns_;
3585      if (modifyCosts&&!badSumPivots) {
3586        int i;
3587        for (i=numberColumns_-1;i>=swapped[iFlip];i--) {
3588          int iSequence=index[i];
3589          double alpha=spare[i];
3590          double value = dj_[iSequence]-theta_*alpha;
3591           
3592          // can't be free here
3593         
3594          if (alpha < 0.0) {
3595            //at upper bound
3596            if (value>dualTolerance_) {
3597              thisIncrease=true;
3598#define MODIFYCOST 2
3599#if MODIFYCOST
3600              // modify cost to hit new tolerance
3601              double modification = alpha*theta_-dj_[iSequence]
3602                +newTolerance;
3603              if ((specialOptions_&(2048+4096+16384))!=0) {
3604                if ((specialOptions_&16384)!=0) {
3605                  if (fabs(modification)<1.0e-8)
3606                    modification=0.0;
3607                } else if ((specialOptions_&2048)!=0) {
3608                  if (fabs(modification)<1.0e-10)
3609                    modification=0.0;
3610                } else {
3611                  if (fabs(modification)<1.0e-12)
3612                    modification=0.0;
3613                }
3614              }
3615              dj_[iSequence] += modification;
3616              cost_[iSequence] +=  modification;
3617              if (modification)
3618                numberChanged_ ++; // Say changed costs
3619              //cost_[iSequence+costOffset] += modification; // save change
3620#endif
3621            }
3622          } else {
3623            // at lower bound
3624            if (-value>dualTolerance_) {
3625              thisIncrease=true;
3626#if MODIFYCOST
3627              // modify cost to hit new tolerance
3628              double modification = alpha*theta_-dj_[iSequence]
3629                -newTolerance;
3630              //modification = CoinMax(modification,-dualTolerance_);
3631              //assert (fabs(modification)<1.0e-7);
3632              if ((specialOptions_&(2048+4096))!=0) {
3633                if ((specialOptions_&2048)!=0) {
3634                  if (fabs(modification)<1.0e-10)
3635                    modification=0.0;
3636                } else {
3637                  if (fabs(modification)<1.0e-12)
3638                    modification=0.0;
3639                }
3640              }
3641              dj_[iSequence] += modification;
3642              cost_[iSequence] +=  modification;
3643              if (modification)
3644                numberChanged_ ++; // Say changed costs
3645              //cost_[iSequence+costOffset] += modification; // save change
3646#endif
3647            }
3648          }
3649        }
3650      }
3651    }
3652  }
3653
3654  if (sequenceIn_>=0) {
3655#ifdef MORE_CAREFUL
3656    // If we have done pivots and things look bad set alpha_ 0.0 to force factorization
3657    if (badSumPivots) {
3658      if (handler_->logLevel()>1)
3659        printf("forcing re-factorization\n");
3660      alpha_=0.0;
3661    }
3662    if (fabs(theta_*badFree)>10.0*dualTolerance_&&factorization_->pivots()) {
3663      if (handler_->logLevel()>1)
3664        printf("forcing re-factorizationon free\n");
3665      alpha_=0.0;
3666    }
3667#endif
3668    lowerIn_ = lower_[sequenceIn_];
3669    upperIn_ = upper_[sequenceIn_];
3670    valueIn_ = solution_[sequenceIn_];
3671    dualIn_ = dj_[sequenceIn_];
3672
3673    if (numberTimesOptimal_) {
3674      // can we adjust cost back closer to original
3675      //*** add coding
3676    }
3677#if MODIFYCOST>1   
3678    // modify cost to hit zero exactly
3679    // so (dualIn_+modification)==theta_*alpha_
3680    double modification = theta_*alpha_-dualIn_;
3681    // But should not move objective too much ??
3682#define DONT_MOVE_OBJECTIVE
3683#ifdef DONT_MOVE_OBJECTIVE
3684    double moveObjective = fabs(modification*solution_[sequenceIn_]);
3685    double smallMove = CoinMax(fabs(objectiveValue_),1.0e-3);
3686    if (moveObjective>smallMove) {
3687      if (handler_->logLevel()>1)
3688        printf("would move objective by %g - original mod %g sol value %g\n",moveObjective,
3689               modification,solution_[sequenceIn_]);
3690      modification *= smallMove/moveObjective;
3691    }
3692#endif
3693    if (badSumPivots)
3694      modification=0.0;
3695    if ((specialOptions_&(2048+4096))!=0) {
3696      if ((specialOptions_&16384)!=0) {
3697        // in fast dual
3698        if (fabs(modification)<1.0e-7)
3699          modification=0.0;
3700      } else if ((specialOptions_&2048)!=0) {
3701        if (fabs(modification)<1.0e-10)
3702          modification=0.0;
3703      } else {
3704        if (fabs(modification)<1.0e-12)
3705          modification=0.0;
3706      }
3707    }
3708    dualIn_ += modification;
3709    dj_[sequenceIn_]=dualIn_;
3710    cost_[sequenceIn_] += modification;
3711    if (modification)
3712      numberChanged_ ++; // Say changed costs
3713    //int costOffset = numberRows_+numberColumns_;
3714    //cost_[sequenceIn_+costOffset] += modification; // save change
3715    //assert (fabs(modification)<1.0e-6);
3716#ifdef CLP_DEBUG
3717    if ((handler_->logLevel()&32)&&fabs(modification)>1.0e-15)
3718      printf("exact %d new cost %g, change %g\n",sequenceIn_,
3719             cost_[sequenceIn_],modification);
3720#endif
3721#endif
3722   
3723    if (alpha_<0.0) {
3724      // as if from upper bound
3725      directionIn_=-1;
3726      upperIn_=valueIn_;
3727    } else {
3728      // as if from lower bound
3729      directionIn_=1;
3730      lowerIn_=valueIn_;
3731    }
3732  }
3733  //if (thisIncrease)
3734  //dualTolerance_+= 1.0e-6*dblParam_[ClpDualTolerance];
3735
3736  // clear arrays
3737
3738  for (i=0;i<2;i++) {
3739    CoinZeroN(array[i],marker[i][0]);
3740    CoinZeroN(array[i]+marker[i][1],numberColumns_-marker[i][1]);
3741  }
3742  return bestPossible;
3743}
3744#ifdef CLP_ALL_ONE_FILE
3745#undef MAXTRY
3746#endif
3747/* Checks if tentative optimal actually means unbounded
3748   Returns -3 if not, 2 if is unbounded */
3749int 
3750ClpSimplexDual::checkUnbounded(CoinIndexedVector * ray,
3751                                   CoinIndexedVector * spare,
3752                                   double changeCost)
3753{
3754  int status=2; // say unbounded
3755  factorization_->updateColumn(spare,ray);
3756  // get reduced cost
3757  int i;
3758  int number=ray->getNumElements();
3759  int * index = ray->getIndices();
3760  double * array = ray->denseVector();
3761  for (i=0;i<number;i++) {
3762    int iRow=index[i];
3763    int iPivot=pivotVariable_[iRow];
3764    changeCost -= cost(iPivot)*array[iRow];
3765  }
3766  double way;
3767  if (changeCost>0.0) {
3768    //try going down
3769    way=1.0;
3770  } else if (changeCost<0.0) {
3771    //try going up
3772    way=-1.0;
3773  } else {
3774#ifdef CLP_DEBUG
3775    printf("can't decide on up or down\n");
3776#endif
3777    way=0.0;
3778    status=-3;
3779  }
3780  double movement=1.0e10*way; // some largish number
3781  double zeroTolerance = 1.0e-14*dualBound_;
3782  for (i=0;i<number;i++) {
3783    int iRow=index[i];
3784    int iPivot=pivotVariable_[iRow];
3785    double arrayValue = array[iRow];
3786    if (fabs(arrayValue)<zeroTolerance)
3787      arrayValue=0.0;
3788    double newValue=solution(iPivot)+movement*arrayValue;
3789    if (newValue>upper(iPivot)+primalTolerance_||
3790        newValue<lower(iPivot)-primalTolerance_)
3791      status=-3; // not unbounded
3792  }
3793  if (status==2) {
3794    // create ray
3795    delete [] ray_;
3796    ray_ = new double [numberColumns_];
3797    CoinZeroN(ray_,numberColumns_);
3798    for (i=0;i<number;i++) {
3799      int iRow=index[i];
3800      int iPivot=pivotVariable_[iRow];
3801      double arrayValue = array[iRow];
3802      if (iPivot<numberColumns_&&fabs(arrayValue)>=zeroTolerance)
3803        ray_[iPivot] = way* array[iRow];
3804    }
3805  }
3806  ray->clear();
3807  return status;
3808}
3809//static int count_alpha=0;
3810/* Checks if finished.  Updates status */
3811void 
3812ClpSimplexDual::statusOfProblemInDual(int & lastCleaned,int type,
3813                                      double * givenDuals, ClpDataSave & saveData,
3814                                      int ifValuesPass)
3815{
3816#ifdef CLP_INVESTIGATE_SERIAL
3817  if (z_thinks>0&&z_thinks<2)
3818    z_thinks+=2;
3819#endif
3820  // If lots of iterations then adjust costs if large ones
3821  if (numberIterations_>4*(numberRows_+numberColumns_)&&objectiveScale_==1.0) {
3822    double largest=0.0;
3823    for (int i=0;i<numberRows_;i++) {
3824      int iColumn = pivotVariable_[i];
3825      largest = CoinMax(largest,fabs(cost_[iColumn]));
3826    }
3827    if (largest>1.0e6) {
3828      objectiveScale_ = 1.0e6/largest;
3829      for (int i=0;i<numberRows_+numberColumns_;i++)
3830        cost_[i] *= objectiveScale_;
3831    }
3832  }
3833  int numberPivots = factorization_->pivots();
3834  double realDualInfeasibilities=0.0;
3835  if (type==2) {
3836    if (alphaAccuracy_!=-1.0)
3837      alphaAccuracy_=-2.0;
3838    // trouble - restore solution
3839    CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
3840    CoinMemcpyN(savedSolution_+numberColumns_ ,
3841           numberRows_,rowActivityWork_);
3842    CoinMemcpyN(savedSolution_ ,
3843           numberColumns_,columnActivityWork_);
3844    // restore extra stuff
3845    int dummy;
3846    matrix_->generalExpanded(this,6,dummy);
3847    forceFactorization_=1; // a bit drastic but ..
3848    changeMade_++; // say something changed
3849    // get correct bounds on all variables
3850    resetFakeBounds(0);
3851  }
3852  int tentativeStatus = problemStatus_;
3853  double changeCost;
3854  bool unflagVariables = true;
3855  bool weightsSaved=false;
3856  bool weightsSaved2=numberIterations_&&!numberPrimalInfeasibilities_;
3857  int dontFactorizePivots = dontFactorizePivots_;
3858  if (type==3) {
3859    type=1;
3860    dontFactorizePivots=1;
3861  }
3862  if (alphaAccuracy_<0.0||!numberPivots||alphaAccuracy_>1.0e4||numberPivots>20) {
3863    if (problemStatus_>-3||numberPivots>dontFactorizePivots) {
3864      // factorize
3865      // later on we will need to recover from singularities
3866      // also we could skip if first time
3867      // save dual weights
3868      dualRowPivot_->saveWeights(this,1);
3869      weightsSaved=true;
3870      if (type) {
3871        // is factorization okay?
3872        if (internalFactorize(1)) {
3873          // no - restore previous basis
3874          unflagVariables = false;
3875          assert (type==1);
3876          changeMade_++; // say something changed
3877          // Keep any flagged variables
3878          int i;
3879          for (i=0;i<numberRows_+numberColumns_;i++) {
3880            if (flagged(i))
3881              saveStatus_[i] |= 64; //say flagged
3882          }
3883          CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
3884          CoinMemcpyN(savedSolution_+numberColumns_ ,
3885                      numberRows_,rowActivityWork_);
3886          CoinMemcpyN(savedSolution_ ,
3887                      numberColumns_,columnActivityWork_);
3888          // restore extra stuff
3889          int dummy;
3890          matrix_->generalExpanded(this,6,dummy);
3891          // get correct bounds on all variables
3892          resetFakeBounds(1);
3893          // need to reject something
3894          char x = isColumn(sequenceOut_) ? 'C' :'R';
3895          handler_->message(CLP_SIMPLEX_FLAG,messages_)
3896            <<x<<sequenceWithin(sequenceOut_)
3897            <<CoinMessageEol;
3898#ifdef COIN_DEVELOP
3899          printf("flag d\n");
3900#endif
3901          setFlagged(sequenceOut_);
3902          progress_.clearBadTimes();
3903         
3904          // Go to safe
3905          factorization_->pivotTolerance(0.99);
3906          forceFactorization_=1; // a bit drastic but ..
3907          type = 2;
3908          //assert (internalFactorize(1)==0);
3909          if (internalFactorize(1)) {
3910            CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
3911            CoinMemcpyN(savedSolution_+numberColumns_ ,
3912                        numberRows_,rowActivityWork_);
3913            CoinMemcpyN(savedSolution_ ,
3914                        numberColumns_,columnActivityWork_);
3915            // restore extra stuff
3916            int dummy;
3917            matrix_->generalExpanded(this,6,dummy);
3918            // debug
3919            int returnCode = internalFactorize(1);
3920            while (returnCode) {
3921              // ouch
3922              // switch off dense
3923              int saveDense = factorization_->denseThreshold();
3924              factorization_->setDenseThreshold(0);
3925              // Go to safe
3926              factorization_->pivotTolerance(0.99);
3927              // make sure will do safe factorization
3928              pivotVariable_[0]=-1;
3929              returnCode=internalFactorize(2);
3930              factorization_->setDenseThreshold(saveDense);
3931            }
3932            // get correct bounds on all variables
3933            resetFakeBounds(1);
3934          }
3935        }
3936      }
3937      if (problemStatus_!=-4||numberPivots>10)
3938        problemStatus_=-3;
3939    }
3940  } else {
3941    //printf("testing with accuracy of %g and status of %d\n",alphaAccuracy_,problemStatus_);
3942    //count_alpha++;
3943    //if ((count_alpha%5000)==0)
3944    //printf("count alpha %d\n",count_alpha);
3945  }
3946  // at this stage status is -3 or -4 if looks infeasible
3947  // get primal and dual solutions
3948#if 0
3949  {
3950    int numberTotal = numberRows_+numberColumns_;
3951    double * saveSol = CoinCopyOfArray(solution_,numberTotal);
3952    double * saveDj = CoinCopyOfArray(dj_,numberTotal);
3953    double tolerance = type ? 1.0e-4 : 1.0e-8;
3954    // always if values pass
3955    double saveObj=objectiveValue_;
3956    double sumPrimal = sumPrimalInfeasibilities_;
3957    int numberPrimal = numberPrimalInfeasibilities_;
3958    double sumDual = sumDualInfeasibilities_;
3959    int numberDual = numberDualInfeasibilities_;
3960    gutsOfSolution(givenDuals,NULL);
3961    int j;
3962    double largestPrimal = tolerance;
3963    int iPrimal=-1;
3964    for (j=0;j<numberTotal;j++) {
3965      double difference = solution_[j]-saveSol[j];
3966      if (fabs(difference)>largestPrimal) {
3967        iPrimal=j;
3968        largestPrimal=fabs(difference);
3969      }
3970    }
3971    double largestDual = tolerance;
3972    int iDual=-1;
3973    for (j=0;j<numberTotal;j++) {
3974      double difference = dj_[j]-saveDj[j];
3975      if (fabs(difference)>largestDual&&upper_[j]>lower_[j]) {
3976        iDual=j;
3977        largestDual=fabs(difference);
3978      }
3979    }
3980    if (!type) {
3981      if (fabs(saveObj-objectiveValue_)>1.0e-5||
3982          numberPrimal!=numberPrimalInfeasibilities_||numberPrimal!=1||
3983          fabs(sumPrimal-sumPrimalInfeasibilities_)>1.0e-5||iPrimal>=0||
3984          numberDual!=numberDualInfeasibilities_||numberDual!=0||
3985          fabs(sumDual-sumDualInfeasibilities_)>1.0e-5||iDual>=0)
3986        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",
3987               type,numberIterations_,numberPivots,
3988               numberPrimal,numberPrimalInfeasibilities_,sumPrimal,sumPrimalInfeasibilities_,
3989               largestPrimal,iPrimal,
3990               numberDual,numberDualInfeasibilities_,sumDual,sumDualInfeasibilities_,
3991               largestDual,iDual,
3992               saveObj,objectiveValue_);
3993    } else {
3994      if (fabs(saveObj-objectiveValue_)>1.0e-5||
3995          numberPrimalInfeasibilities_||iPrimal>=0||
3996          numberDualInfeasibilities_||iDual>=0)
3997        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",
3998               type,numberIterations_,numberPivots,
3999               numberPrimal,numberPrimalInfeasibilities_,sumPrimal,sumPrimalInfeasibilities_,
4000               largestPrimal,iPrimal,
4001               numberDual,numberDualInfeasibilities_,sumDual,sumDualInfeasibilities_,
4002               largestDual,iDual,
4003               saveObj,objectiveValue_);
4004    }
4005    delete [] saveSol;
4006    delete [] saveDj;
4007  }
4008#else
4009  if (type||ifValuesPass)
4010    gutsOfSolution(givenDuals,NULL);
4011#endif
4012  // If bad accuracy treat as singular
4013  if ((largestPrimalError_>1.0e15||largestDualError_>1.0e15)&&numberIterations_) {
4014    // restore previous basis
4015    unflagVariables = false;
4016    changeMade_++; // say something changed
4017    // Keep any flagged variables
4018    int i;
4019    for (i=0;i<numberRows_+numberColumns_;i++) {
4020      if (flagged(i))
4021        saveStatus_[i] |= 64; //say flagged
4022    }
4023    CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
4024    CoinMemcpyN(savedSolution_+numberColumns_ ,
4025           numberRows_,rowActivityWork_);
4026    CoinMemcpyN(savedSolution_ ,
4027           numberColumns_,columnActivityWork_);
4028    // restore extra stuff
4029    int dummy;
4030    matrix_->generalExpanded(this,6,dummy);
4031    // get correct bounds on all variables
4032    resetFakeBounds(1);
4033    // need to reject something
4034    char x = isColumn(sequenceOut_) ? 'C' :'R';
4035    handler_->message(CLP_SIMPLEX_FLAG,messages_)
4036      <<x<<sequenceWithin(sequenceOut_)
4037      <<CoinMessageEol;
4038#ifdef COIN_DEVELOP
4039    printf("flag e\n");
4040#endif
4041    setFlagged(sequenceOut_);
4042    progress_.clearBadTimes();
4043   
4044    // Go to safer
4045    double newTolerance = CoinMin(1.1*factorization_->pivotTolerance(),0.99);
4046    factorization_->pivotTolerance(newTolerance);
4047    forceFactorization_=1; // a bit drastic but ..
4048    if (alphaAccuracy_!=-1.0)
4049      alphaAccuracy_=-2.0;
4050    type = 2;
4051    //assert (internalFactorize(1)==0);
4052    if (internalFactorize(1)) {
4053      CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
4054      CoinMemcpyN(savedSolution_+numberColumns_ ,
4055              numberRows_,rowActivityWork_);
4056      CoinMemcpyN(savedSolution_ ,
4057           numberColumns_,columnActivityWork_);
4058      // restore extra stuff
4059      int dummy;
4060      matrix_->generalExpanded(this,6,dummy);
4061      // debug
4062      int returnCode = internalFactorize(1);
4063      while (returnCode) {
4064        // ouch
4065        // switch off dense
4066        int saveDense = factorization_->denseThreshold();
4067        factorization_->setDenseThreshold(0);
4068        // Go to safe
4069        factorization_->pivotTolerance(0.99);
4070        // make sure will do safe factorization
4071        pivotVariable_[0]=-1;
4072        returnCode=internalFactorize(2);
4073        factorization_->setDenseThreshold(saveDense);
4074      }
4075      // get correct bounds on all variables
4076      resetFakeBounds(1);
4077    }
4078    // get primal and dual solutions
4079    gutsOfSolution(givenDuals,NULL);
4080  } else if (goodAccuracy()) {
4081    // Can reduce tolerance
4082    double newTolerance = CoinMax(0.99*factorization_->pivotTolerance(),saveData.pivotTolerance_);
4083    factorization_->pivotTolerance(newTolerance);
4084  }
4085  bool reallyBadProblems=false;
4086  // Double check infeasibility if no action
4087  if (progress_.lastIterationNumber(0)==numberIterations_) {
4088    if (dualRowPivot_->looksOptimal()) {
4089      numberPrimalInfeasibilities_ = 0;
4090      sumPrimalInfeasibilities_ = 0.0;
4091    }
4092#if 1
4093  } else {
4094    double thisObj = objectiveValue_-bestPossibleImprovement_;
4095#ifdef CLP_INVESTIGATE
4096    assert (bestPossibleImprovement_>-1000.0&&objectiveValue_>-1.0e100);
4097    if (bestPossibleImprovement_)
4098      printf("obj %g add in %g -> %g\n",objectiveValue_,bestPossibleImprovement_,
4099             thisObj);
4100#endif
4101    double lastObj = progress_.lastObjective(0);
4102#ifndef NDEBUG
4103#ifdef COIN_DEVELOP
4104    resetFakeBounds(-1);
4105#endif
4106#endif
4107#ifdef CLP_REPORT_PROGRESS
4108    ixxxxxx++;
4109    if (ixxxxxx>=ixxyyyy-4&&ixxxxxx<=ixxyyyy) {
4110      char temp[20];
4111      sprintf(temp,"sol%d.out",ixxxxxx);
4112      printf("sol%d.out\n",ixxxxxx);
4113      FILE * fp = fopen(temp,"w");
4114      int nTotal=numberRows_+numberColumns_;
4115      for (int i=0;i<nTotal;i++)
4116        fprintf(fp,"%d %d %g %g %g %g %g\n",
4117                i,status_[i],lower_[i],solution_[i],upper_[i],cost_[i],dj_[i]);
4118      fclose(fp);
4119    }
4120#endif
4121    if(!ifValuesPass&&firstFree_<0) {
4122      double testTol = 5.0e-3;
4123      if (progress_.timesFlagged()>10) {
4124        testTol *= pow(2.0,progress_.timesFlagged()-8);
4125      } else if (progress_.timesFlagged()>5) {
4126        testTol *=5.0;
4127      }
4128      if (lastObj>thisObj+
4129          testTol*(fabs(thisObj)+fabs(lastObj))+testTol) {
4130        int maxFactor = factorization_->maximumPivots();
4131        if ((specialOptions_&1048576)==0) {
4132          if (progress_.timesFlagged()>10) 
4133            progress_.incrementReallyBadTimes();
4134          if (maxFactor>10-9) {
4135#ifdef COIN_DEVELOP
4136            printf("lastobj %g thisobj %g\n",lastObj,thisObj);
4137#endif
4138            //if (forceFactorization_<0)
4139            //forceFactorization_= maxFactor;
4140            //forceFactorization_ = CoinMax(1,(forceFactorization_>>1));
4141            if ((progressFlag_&4)==0&&lastObj<thisObj+1.0e4&&
4142                largestPrimalError_<1.0e2) {
4143              // Just save costs
4144              // save extra copy of cost_
4145              int nTotal=numberRows_+numberColumns_;
4146              double * temp = new double [2*nTotal];
4147              memcpy(temp,cost_,nTotal*sizeof(double));
4148              memcpy(temp+nTotal,cost_,nTotal*sizeof(double));
4149              delete [] cost_;
4150              cost_=temp;
4151              objectiveWork_ = cost_;
4152              rowObjectiveWork_ = cost_+numberColumns_;
4153              progressFlag_ |= 4;
4154            } else {
4155              forceFactorization_=1;
4156#ifdef COIN_DEVELOP
4157              printf("Reducing factorization frequency - bad backwards\n");
4158#endif
4159#if 1
4160              unflagVariables = false;
4161              changeMade_++; // say something changed
4162              int nTotal = numberRows_+numberColumns_;
4163              CoinMemcpyN(saveStatus_,nTotal,status_);
4164              CoinMemcpyN(savedSolution_+numberColumns_ ,
4165                          numberRows_,rowActivityWork_);
4166              CoinMemcpyN(savedSolution_ ,
4167                          numberColumns_,columnActivityWork_);
4168              if ((progressFlag_&4)==0) {
4169                // save extra copy of cost_
4170                double * temp = new double [2*nTotal];
4171                memcpy(temp,cost_,nTotal*sizeof(double));
4172                memcpy(temp+nTotal,cost_,nTotal*sizeof(double));
4173                delete [] cost_;
4174                cost_=temp;
4175                objectiveWork_ = cost_;
4176                rowObjectiveWork_ = cost_+numberColumns_;
4177                progressFlag_ |= 4;
4178              } else {
4179                memcpy(cost_,cost_+nTotal,nTotal*sizeof(double));
4180              }
4181              // restore extra stuff
4182              int dummy;
4183              matrix_->generalExpanded(this,6,dummy);
4184              double pivotTolerance = factorization_->pivotTolerance();
4185              if(pivotTolerance<0.2)
4186                factorization_->pivotTolerance(0.2);
4187              else if(progress_.timesFlagged()>2)
4188                factorization_->pivotTolerance(CoinMin(pivotTolerance*1.1,0.99));
4189              if (alphaAccuracy_!=-1.0)
4190                alphaAccuracy_=-2.0;
4191              if (internalFactorize(1)) {
4192                CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
4193                CoinMemcpyN(savedSolution_+numberColumns_ ,
4194                            numberRows_,rowActivityWork_);
4195                CoinMemcpyN(savedSolution_ ,
4196                            numberColumns_,columnActivityWork_);
4197                // restore extra stuff
4198                int dummy;
4199                matrix_->generalExpanded(this,6,dummy);
4200                // debug
4201                int returnCode = internalFactorize(1);
4202                while (returnCode) {
4203                  // ouch
4204                  // switch off dense
4205                  int saveDense = factorization_->denseThreshold();
4206                  factorization_->setDenseThreshold(0);
4207                  // Go to safe
4208                  factorization_->pivotTolerance(0.99);
4209                  // make sure will do safe factorization
4210                  pivotVariable_[0]=-1;
4211                  returnCode=internalFactorize(2);
4212                  factorization_->setDenseThreshold(saveDense);
4213                }
4214              }
4215              resetFakeBounds(0);
4216              type = 2; // so will restore weights
4217              // get primal and dual solutions
4218              gutsOfSolution(givenDuals,NULL);
4219              if (numberPivots<2) {
4220                // need to reject something
4221                char x = isColumn(sequenceOut_) ? 'C' :'R';
4222                handler_->message(CLP_SIMPLEX_FLAG,messages_)
4223                  <<x<<sequenceWithin(sequenceOut_)
4224                  <<CoinMessageEol;
4225#ifdef COIN_DEVELOP
4226                printf("flag d\n");
4227#endif
4228                setFlagged(sequenceOut_);
4229                progress_.clearBadTimes();
4230                progress_.incrementTimesFlagged();
4231              }
4232              if (numberPivots<10)
4233                reallyBadProblems=true;
4234#ifdef COIN_DEVELOP
4235              printf("obj now %g\n",objectiveValue_);
4236#endif
4237              progress_.modifyObjective(objectiveValue_
4238                                        -bestPossibleImprovement_);
4239#endif
4240            }
4241          }
4242        } else {
4243          // in fast dual give up
4244#ifdef COIN_DEVELOP
4245          printf("In fast dual?\n");
4246#endif
4247          problemStatus_=3;
4248        }
4249      } else if (lastObj<thisObj-1.0e-5*CoinMax(fabs(thisObj),fabs(lastObj))-1.0e-3) {
4250        numberTimesOptimal_=0;
4251      }
4252    }
4253#endif
4254  }
4255  // Up tolerance if looks a bit odd
4256  if (numberIterations_>CoinMax(1000,numberRows_>>4)&&(specialOptions_&64)!=0) {
4257    if (sumPrimalInfeasibilities_&&sumPrimalInfeasibilities_<1.0e5) {
4258      int backIteration = progress_.lastIterationNumber(CLP_PROGRESS-1);
4259      if (backIteration>0&&numberIterations_-backIteration<9*CLP_PROGRESS) {
4260        if (factorization_->pivotTolerance()<0.9) {
4261          // up tolerance
4262          factorization_->pivotTolerance(CoinMin(factorization_->pivotTolerance()*1.05+0.02,0.91));
4263          //printf("tol now %g\n",factorization_->pivotTolerance());
4264          progress_.clearIterationNumbers();
4265        }
4266      }
4267    }
4268  }
4269  // Check if looping
4270  int loop;
4271  if (!givenDuals&&type!=2) 
4272    loop = progress_.looping();
4273  else
4274    loop=-1;
4275  if (progress_.reallyBadTimes()>10) {
4276    problemStatus_ = 10; // instead - try other algorithm
4277#if COIN_DEVELOP>2
4278    printf("returning at %d\n",__LINE__);
4279#endif
4280  }
4281  int situationChanged=0;
4282  if (loop>=0) {
4283    problemStatus_ = loop; //exit if in loop
4284    if (!problemStatus_) {
4285      // declaring victory
4286      numberPrimalInfeasibilities_ = 0;
4287      sumPrimalInfeasibilities_ = 0.0;
4288    } else {
4289      problemStatus_ = 10; // instead - try other algorithm
4290#if COIN_DEVELOP>2
4291      printf("returning at %d\n",__LINE__);
4292#endif
4293    }
4294    return;
4295  } else if (loop<-1) {
4296    // something may have changed
4297    gutsOfSolution(NULL,NULL);
4298    situationChanged=1;
4299  }
4300  // really for free variables in
4301  if((progressFlag_&2)!=0) {
4302    situationChanged=2;
4303  }
4304  progressFlag_ &= (~3); //reset progress flag
4305  if ((progressFlag_&4)!=0) {
4306    // save copy of cost_
4307    int nTotal = numberRows_+numberColumns_;
4308    memcpy(cost_+nTotal,cost_,nTotal*sizeof(double));
4309  }
4310  /*if (!numberIterations_&&sumDualInfeasibilities_)
4311    printf("OBJ %g sumPinf %g sumDinf %g\n",
4312           objectiveValue(),sumPrimalInfeasibilities_,
4313           sumDualInfeasibilities_);*/
4314  if (handler_->detail(CLP_SIMPLEX_STATUS,messages_)<100) {
4315    handler_->message(CLP_SIMPLEX_STATUS,messages_)
4316      <<numberIterations_<<objectiveValue();
4317    handler_->printing(sumPrimalInfeasibilities_>0.0)
4318      <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
4319    handler_->printing(sumDualInfeasibilities_>0.0)
4320      <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
4321    handler_->printing(numberDualInfeasibilitiesWithoutFree_
4322                       <numberDualInfeasibilities_)
4323                         <<numberDualInfeasibilitiesWithoutFree_;
4324    handler_->message()<<CoinMessageEol;
4325  }
4326#ifdef CLP_REPORT_PROGRESS
4327    if (ixxxxxx>=ixxyyyy-4&&ixxxxxx<=ixxyyyy) {
4328      char temp[20];
4329      sprintf(temp,"x_sol%d.out",ixxxxxx);
4330      FILE * fp = fopen(temp,"w");
4331      int nTotal=numberRows_+numberColumns_;
4332      for (int i=0;i<nTotal;i++)
4333        fprintf(fp,"%d %d %g %g %g %g %g\n",
4334                i,status_[i],lower_[i],solution_[i],upper_[i],cost_[i],dj_[i]);
4335      fclose(fp);
4336      if (ixxxxxx==ixxyyyy)
4337        exit(6);
4338    }
4339#endif
4340  realDualInfeasibilities=sumDualInfeasibilities_;
4341  double saveTolerance =dualTolerance_;
4342  // If we need to carry on cleaning variables
4343  if (!numberPrimalInfeasibilities_&&(specialOptions_&1024)!=0&&CLEAN_FIXED) {
4344    for (int iRow=0;iRow<numberRows_;iRow++) {
4345      int iPivot=pivotVariable_[iRow];
4346      if (!flagged(iPivot)&&pivoted(iPivot)) {
4347        // carry on
4348        numberPrimalInfeasibilities_=-1;
4349        sumOfRelaxedPrimalInfeasibilities_ = 1.0;
4350        sumPrimalInfeasibilities_ = 1.0;
4351        break;
4352      }
4353    }
4354  }
4355  /* If we are primal feasible and any dual infeasibilities are on
4356     free variables then it is better to go to primal */
4357  if (!numberPrimalInfeasibilities_&&!numberDualInfeasibilitiesWithoutFree_&&
4358      numberDualInfeasibilities_)
4359    problemStatus_=10;
4360  // dual bound coming in
4361  //double saveDualBound = dualBound_;
4362  bool needCleanFake=false;
4363  while (problemStatus_<=-3) {
4364    int cleanDuals=0;
4365    if (situationChanged!=0)
4366      cleanDuals=1;
4367    int numberChangedBounds=0;
4368    int doOriginalTolerance=0;
4369    if ( lastCleaned==numberIterations_)
4370      doOriginalTolerance=1;
4371    // check optimal
4372    // give code benefit of doubt
4373    if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
4374        sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
4375      // say optimal (with these bounds etc)
4376      numberDualInfeasibilities_ = 0;
4377      sumDualInfeasibilities_ = 0.0;
4378      numberPrimalInfeasibilities_ = 0;
4379      sumPrimalInfeasibilities_ = 0.0;
4380    }
4381    //if (dualFeasible()||problemStatus_==-4||(primalFeasible()&&!numberDualInfeasibilitiesWithoutFree_)) {
4382    if (dualFeasible()||problemStatus_==-4) {
4383      progress_.modifyObjective(objectiveValue_
4384                                -bestPossibleImprovement_);
4385#ifdef COIN_DEVELOP
4386      if (sumDualInfeasibilities_||bestPossibleImprovement_)
4387        printf("improve %g dualinf %g -> %g\n",
4388               bestPossibleImprovement_,sumDualInfeasibilities_,
4389               sumDualInfeasibilities_*dualBound_);
4390#endif
4391      // see if cutoff reached
4392      double limit = 0.0;
4393      getDblParam(ClpDualObjectiveLimit, limit);
4394#if 0
4395      if(fabs(limit)<1.0e30&&objectiveValue()*optimizationDirection_>
4396         limit+1.0e-7+1.0e-8*fabs(limit)&&!numberAtFakeBound()) {
4397        //looks infeasible on objective
4398        if (perturbation_==101) {
4399          cleanDuals=1;
4400          // Save costs
4401          int numberTotal = numberRows_+numberColumns_;
4402          double * saveCost = CoinCopyOfArray(cost_,numberTotal);
4403          // make sure fake bounds are back
4404          changeBounds(1,NULL,changeCost);
4405          createRim4(false);
4406          // make sure duals are current
4407          computeDuals(givenDuals);
4408          checkDualSolution();
4409          if(objectiveValue()*optimizationDirection_>
4410             limit+1.0e-7+1.0e-8*fabs(limit)&&!numberDualInfeasibilities_) {
4411            perturbation_=102; // stop any perturbations
4412            printf("cutoff test succeeded\n");
4413          } else {
4414            printf("cutoff test failed\n");
4415            // put back
4416            memcpy(cost_,saveCost,numberTotal*sizeof(double));
4417            // make sure duals are current
4418            computeDuals(givenDuals);
4419            checkDualSolution();
4420            progress_.modifyObjective(-COIN_DBL_MAX);
4421            problemStatus_=-1;
4422          }
4423          delete [] saveCost;
4424        }
4425      }
4426#endif
4427      if (primalFeasible()&&!givenDuals) {
4428        // may be optimal - or may be bounds are wrong
4429        handler_->message(CLP_DUAL_BOUNDS,messages_)
4430          <<dualBound_
4431          <<CoinMessageEol;
4432        // save solution in case unbounded
4433        double * saveColumnSolution = NULL;
4434        double * saveRowSolution = NULL;
4435        bool inCbc = (specialOptions_&(0x01000000|16384))!=0;
4436        if (!inCbc) {
4437          saveColumnSolution = CoinCopyOfArray(columnActivityWork_,numberColumns_);
4438          saveRowSolution = CoinCopyOfArray(rowActivityWork_,numberRows_);
4439        }
4440        numberChangedBounds=changeBounds(0,rowArray_[3],changeCost);
4441        if (numberChangedBounds<=0&&!numberDualInfeasibilities_) {
4442          //looks optimal - do we need to reset tolerance
4443          if (perturbation_==101) {
4444            perturbation_=102; // stop any perturbations
4445            cleanDuals=1;
4446            // make sure fake bounds are back
4447            //computeObjectiveValue();
4448            changeBounds(1,NULL,changeCost);
4449            //computeObjectiveValue();
4450            createRim4(false);
4451            // make sure duals are current
4452            computeDuals(givenDuals);
4453            checkDualSolution();
4454            progress_.modifyObjective(-COIN_DBL_MAX);
4455#define DUAL_TRY_FASTER
4456#ifdef DUAL_TRY_FASTER
4457            if (numberDualInfeasibilities_) {
4458#endif
4459              numberChanged_=1; // force something to happen
4460              lastCleaned=numberIterations_-1;
4461#ifdef DUAL_TRY_FASTER
4462            } else {
4463              //double value = objectiveValue_;
4464              computeObjectiveValue(true);
4465              //printf("old %g new %g\n",value,objectiveValue_);
4466              //numberChanged_=1;
4467            }
4468#endif
4469          }
4470          if (lastCleaned<numberIterations_&&numberTimesOptimal_<4&&
4471              (numberChanged_||(specialOptions_&4096)==0)) {
4472            doOriginalTolerance=2;
4473            numberTimesOptimal_++;
4474            changeMade_++; // say something changed
4475            if (numberTimesOptimal_==1) {
4476              dualTolerance_ = dblParam_[ClpDualTolerance];
4477            } else {
4478              if (numberTimesOptimal_==2) {
4479                // better to have small tolerance even if slower
4480                factorization_->zeroTolerance(CoinMin(factorization_->zeroTolerance(),1.0e-15));
4481              }
4482              dualTolerance_ = dblParam_[ClpDualTolerance];
4483              dualTolerance_ *= pow(2.0,numberTimesOptimal_-1);
4484            }
4485            cleanDuals=2; // If nothing changed optimal else primal
4486          } else {
4487            problemStatus_=0; // optimal
4488            if (lastCleaned<numberIterations_&&numberChanged_) {
4489              handler_->message(CLP_SIMPLEX_GIVINGUP,messages_)
4490                <<CoinMessageEol;
4491            }
4492          }
4493        } else {
4494          cleanDuals=1;
4495          if (doOriginalTolerance==1) {
4496            // check unbounded
4497            // find a variable with bad dj
4498            int iSequence;
4499            int iChosen=-1;
4500            if (!inCbc) {
4501              double largest = 100.0*primalTolerance_;
4502              for (iSequence=0;iSequence<numberRows_+numberColumns_;
4503                   iSequence++) {
4504                double djValue = dj_[iSequence];
4505                double originalLo = originalLower(iSequence);
4506                double originalUp = originalUpper(iSequence);
4507                if (fabs(djValue)>fabs(largest)) {
4508                  if (getStatus(iSequence)!=basic) {
4509                    if (djValue>0&&originalLo<-1.0e20) {
4510                      if (djValue>fabs(largest)) {
4511                        largest=djValue;
4512                        iChosen=iSequence;
4513                      }
4514                    } else if (djValue<0&&originalUp>1.0e20) {
4515                      if (-djValue>fabs(largest)) {
4516                        largest=djValue;
4517                        iChosen=iSequence;
4518                      }
4519                    } 
4520                  }
4521                }
4522              }
4523            }
4524            if (iChosen>=0) {
4525              int iSave=sequenceIn_;
4526              sequenceIn_=iChosen;
4527              unpack(rowArray_[1]);
4528              sequenceIn_ = iSave;
4529              // if dual infeasibilities then must be free vector so add in dual
4530              if (numberDualInfeasibilities_) {
4531                if (fabs(changeCost)>1.0e-5)
4532                  printf("Odd free/unbounded combo\n");
4533                changeCost += cost_[iChosen];
4534              }
4535              problemStatus_ = checkUnbounded(rowArray_[1],rowArray_[0],
4536                                              changeCost);
4537              rowArray_[1]->clear();
4538            } else {
4539              problemStatus_=-3;
4540            }
4541            if (problemStatus_==2&&perturbation_==101) {
4542              perturbation_=102; // stop any perturbations
4543              cleanDuals=1;
4544              createRim4(false);
4545              progress_.modifyObjective(-COIN_DBL_MAX);
4546              problemStatus_=-1;
4547            }
4548            if (problemStatus_==2) {
4549              // it is unbounded - restore solution
4550              // but first add in changes to non-basic
4551              int iColumn;
4552              double * original = columnArray_[0]->denseVector();
4553              for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4554                if(getColumnStatus(iColumn)!= basic)
4555                  ray_[iColumn] += 
4556                    saveColumnSolution[iColumn]-original[iColumn];
4557                columnActivityWork_[iColumn] = original[iColumn];
4558              }
4559              CoinMemcpyN(saveRowSolution,numberRows_,
4560                                rowActivityWork_);
4561            }
4562          } else {
4563            doOriginalTolerance=2;
4564            rowArray_[0]->clear();
4565          }
4566        }
4567        delete [] saveColumnSolution;
4568        delete [] saveRowSolution;
4569      }
4570      if (problemStatus_==-4||problemStatus_==-5) {
4571        // may be infeasible - or may be bounds are wrong
4572        numberChangedBounds=changeBounds(0,NULL,changeCost);
4573        needCleanFake=true;
4574        /* Should this be here as makes no difference to being feasible.
4575           But seems to make a difference to run times. */
4576        if (perturbation_==101&&0) {
4577          perturbation_=102; // stop any perturbations
4578          cleanDuals=1;
4579          numberChangedBounds=1;
4580          // make sure fake bounds are back
4581          changeBounds(1,NULL,changeCost);
4582          needCleanFake=true;
4583          createRim4(false);
4584          progress_.modifyObjective(-COIN_DBL_MAX);
4585        }
4586        if ((numberChangedBounds<=0||dualBound_>1.0e20||
4587            (largestPrimalError_>1.0&&dualBound_>1.0e17))&&
4588            (numberPivots<4||sumPrimalInfeasibilities_>1.0e-6)) {
4589          problemStatus_=1; // infeasible
4590          if (perturbation_==101) {
4591            perturbation_=102; // stop any perturbations
4592            //cleanDuals=1;
4593            //numberChangedBounds=1;
4594            //createRim4(false);
4595          }
4596        } else {
4597          problemStatus_=-1; //iterate
4598          cleanDuals=1;
4599          if (numberChangedBounds<=0)
4600            doOriginalTolerance=2;
4601          // and delete ray which has been created
4602          delete [] ray_;
4603          ray_ = NULL;
4604        }
4605
4606      }
4607    } else {
4608      cleanDuals=1;
4609    }
4610    if (problemStatus_<0) {
4611      if (doOriginalTolerance==2) {
4612        // put back original tolerance
4613        lastCleaned=numberIterations_;
4614        numberChanged_ =0; // Number of variables with changed costs
4615        handler_->message(CLP_DUAL_ORIGINAL,messages_)
4616          <<CoinMessageEol;
4617        perturbation_=102; // stop any perturbations
4618#if 0
4619        double * xcost = new double[numberRows_+numberColumns_];
4620        double * xlower = new double[numberRows_+numberColumns_];
4621        double * xupper = new double[numberRows_+numberColumns_];
4622        double * xdj = new double[numberRows_+numberColumns_];
4623        double * xsolution = new double[numberRows_+numberColumns_];
4624 CoinMemcpyN(cost_,(numberRows_+numberColumns_),xcost);
4625 CoinMemcpyN(lower_,(numberRows_+numberColumns_),xlower);
4626 CoinMemcpyN(upper_,(numberRows_+numberColumns_),xupper);
4627 CoinMemcpyN(dj_,(numberRows_+numberColumns_),xdj);
4628 CoinMemcpyN(solution_,(numberRows_+numberColumns_),xsolution);
4629#endif
4630        createRim4(false);
4631        progress_.modifyObjective(-COIN_DBL_MAX);
4632        // make sure duals are current
4633        computeDuals(givenDuals);
4634        checkDualSolution();
4635#if 0
4636        int i;
4637        for (i=0;i<numberRows_+numberColumns_;i++) {
4638          if (cost_[i]!=xcost[i])
4639            printf("** %d old cost %g new %g sol %g\n",
4640                   i,xcost[i],cost_[i],solution_[i]);
4641          if (lower_[i]!=xlower[i])
4642            printf("** %d old lower %g new %g sol %g\n",
4643                   i,xlower[i],lower_[i],solution_[i]);
4644          if (upper_[i]!=xupper[i])
4645            printf("** %d old upper %g new %g sol %g\n",
4646                   i,xupper[i],upper_[i],solution_[i]);
4647          if (dj_[i]!=xdj[i])
4648            printf("** %d old dj %g new %g sol %g\n",
4649                   i,xdj[i],dj_[i],solution_[i]);
4650          if (solution_[i]!=xsolution[i])
4651            printf("** %d old solution %g new %g sol %g\n",
4652                   i,xsolution[i],solution_[i],solution_[i]);
4653        }
4654        //delete [] xcost;
4655        //delete [] xupper;
4656        //delete [] xlower;
4657        //delete [] xdj;
4658        //delete [] xsolution;
4659#endif
4660        // put back bounds as they were if was optimal
4661        if (doOriginalTolerance==2&&cleanDuals!=2) {
4662          changeMade_++; // say something changed
4663          /* We may have already changed some bounds in this function
4664             so save numberFake_ and add in.
4665
4666             Worst that can happen is that we waste a bit of time  - but it must be finite.
4667          */
4668          //int saveNumberFake = numberFake_;
4669          //resetFakeBounds(-1);
4670          changeBounds(3,NULL,changeCost);
4671          needCleanFake=true;
4672          //numberFake_ += saveNumberFake;
4673          //resetFakeBounds(-1);
4674          cleanDuals=2;
4675          //cleanDuals=1;
4676        }
4677#if 0
4678        //int i;
4679        for (i=0;i<numberRows_+numberColumns_;i++) {
4680          if (cost_[i]!=xcost[i])
4681            printf("** %d old cost %g new %g sol %g\n",
4682                   i,xcost[i],cost_[i],solution_[i]);
4683          if (lower_[i]!=xlower[i])
4684            printf("** %d old lower %g new %g sol %g\n",
4685                   i,xlower[i],lower_[i],solution_[i]);
4686          if (upper_[i]!=xupper[i])
4687            printf("** %d old upper %g new %g sol %g\n",
4688                   i,xupper[i],upper_[i],solution_[i]);
4689          if (dj_[i]!=xdj[i])
4690            printf("** %d old dj %g new %g sol %g\n",
4691                   i,xdj[i],dj_[i],solution_[i]);
4692          if (solution_[i]!=xsolution[i])
4693            printf("** %d old solution %g new %g sol %g\n",
4694                   i,xsolution[i],solution_[i],solution_[i]);
4695        }
4696        delete [] xcost;
4697        delete [] xupper;
4698        delete [] xlower;
4699        delete [] xdj;
4700        delete [] xsolution;
4701#endif
4702      }
4703      if (cleanDuals==1||(cleanDuals==2&&!numberDualInfeasibilities_)) {
4704        // make sure dual feasible
4705        // look at all rows and columns
4706        rowArray_[0]->clear();
4707        columnArray_[0]->clear();
4708        double objectiveChange=0.0;
4709#if 0
4710        double * xcost = new double[numberRows_+numberColumns_];
4711        double * xlower = new double[numberRows_+numberColumns_];
4712        double * xupper = new double[numberRows_+numberColumns_];
4713        double * xdj = new double[numberRows_+numberColumns_];
4714        double * xsolution = new double[numberRows_+numberColumns_];
4715 CoinMemcpyN(cost_,(numberRows_+numberColumns_),xcost);
4716 CoinMemcpyN(lower_,(numberRows_+numberColumns_),xlower);
4717 CoinMemcpyN(upper_,(numberRows_+numberColumns_),xupper);
4718 CoinMemcpyN(dj_,(numberRows_+numberColumns_),xdj);
4719 CoinMemcpyN(solution_,(numberRows_+numberColumns_),xsolution);
4720#endif
4721        if (givenDuals)
4722          dualTolerance_=1.0e50;
4723        updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[1],
4724          0.0,objectiveChange,true);
4725        dualTolerance_=saveTolerance;
4726#if 0
4727        int i;
4728        for (i=0;i<numberRows_+numberColumns_;i++) {
4729          if (cost_[i]!=xcost[i])
4730            printf("** %d old cost %g new %g sol %g\n",
4731                   i,xcost[i],cost_[i],solution_[i]);
4732          if (lower_[i]!=xlower[i])
4733            printf("** %d old lower %g new %g sol %g\n",
4734                   i,xlower[i],lower_[i],solution_[i]);
4735          if (upper_[i]!=xupper[i])
4736            printf("** %d old upper %g new %g sol %g\n",
4737                   i,xupper[i],upper_[i],solution_[i]);
4738          if (dj_[i]!=xdj[i])
4739            printf("** %d old dj %g new %g sol %g\n",
4740                   i,xdj[i],dj_[i],solution_[i]);
4741          if (solution_[i]!=xsolution[i])
4742            printf("** %d old solution %g new %g sol %g\n",
4743                   i,xsolution[i],solution_[i],solution_[i]);
4744        }
4745        delete [] xcost;
4746        delete [] xupper;
4747        delete [] xlower;
4748        delete [] xdj;
4749        delete [] xsolution;
4750#endif
4751        // for now - recompute all
4752        gutsOfSolution(NULL,NULL);
4753        if (givenDuals)
4754          dualTolerance_=1.0e50;
4755        updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[1],
4756          0.0,objectiveChange,true);
4757        dualTolerance_=saveTolerance;
4758        //assert(numberDualInfeasibilitiesWithoutFree_==0);
4759        if (numberDualInfeasibilities_) {
4760          if (numberPrimalInfeasibilities_||numberPivots) {
4761            problemStatus_=-1; // carry on as normal
4762          } else {
4763            problemStatus_=10; // try primal
4764#if COIN_DEVELOP>1
4765            printf("returning at %d\n",__LINE__);
4766#endif
4767          }
4768        } else if (situationChanged==2) {
4769          problemStatus_=-1; // carry on as normal
4770        }
4771        situationChanged=0;
4772      } else {
4773        // iterate
4774        if (cleanDuals!=2) {
4775          problemStatus_=-1;
4776        } else {
4777          problemStatus_ = 10; // try primal
4778#if COIN_DEVELOP>2
4779          printf("returning at %d\n",__LINE__);
4780#endif
4781        }
4782      }
4783    }
4784  }
4785  // unflag all variables (we may want to wait a bit?)
4786  if ((tentativeStatus!=-2&&tentativeStatus!=-1)&&unflagVariables) {
4787    int iRow;
4788    int numberFlagged=0;
4789    for (iRow=0;iRow<numberRows_;iRow++) {
4790      int iPivot=pivotVariable_[iRow];
4791      if (flagged(iPivot)) {
4792        numberFlagged++;
4793        clearFlagged(iPivot);
4794      }
4795    }
4796#ifdef COIN_DEVELOP
4797    if (numberFlagged) {
4798      printf("unflagging %d variables - tentativeStatus %d probStat %d ninf %d nopt %d\n",numberFlagged,tentativeStatus,
4799             problemStatus_,numberPrimalInfeasibilities_,
4800             numberTimesOptimal_);
4801    }
4802#endif
4803    unflagVariables = numberFlagged>0;
4804    if (numberFlagged&&!numberPivots) {
4805      /* looks like trouble as we have not done any iterations.
4806         Try changing pivot tolerance then give it a few goes and give up */
4807      if (factorization_->pivotTolerance()<0.9) {
4808        factorization_->pivotTolerance(0.99);
4809        problemStatus_=-1;
4810      } else if (numberTimesOptimal_<3) {
4811        numberTimesOptimal_++;
4812        problemStatus_=-1;
4813      } else {
4814        unflagVariables=false;
4815        //secondaryStatus_ = 1; // and say probably infeasible
4816        // try primal
4817        problemStatus_=10;
4818#if COIN_DEVELOP>1
4819        printf("returning at %d\n",__LINE__);
4820#endif
4821      }
4822    }
4823  }
4824  if (problemStatus_<0) {
4825    if (needCleanFake) {
4826      double dummyChangeCost=0.0;
4827      changeBounds(3,NULL,dummyChangeCost);
4828    }
4829#if 0
4830    if (objectiveValue_<lastObjectiveValue_-1.0e-8*
4831        CoinMax(fabs(objectivevalue_),fabs(lastObjectiveValue_))) {
4832    } else {
4833      lastObjectiveValue_=objectiveValue_;
4834    }
4835#endif
4836    if (type==0||type==1) {
4837      if (!type) {
4838        // create save arrays
4839        delete [] saveStatus_;
4840        delete [] savedSolution_;
4841        saveStatus_ = new unsigned char [numberRows_+numberColumns_];
4842        savedSolution_ = new double [numberRows_+numberColumns_];
4843      }
4844      // save arrays
4845      CoinMemcpyN(status_,numberColumns_+numberRows_,saveStatus_);
4846      CoinMemcpyN(rowActivityWork_,
4847                  numberRows_,savedSolution_+numberColumns_);
4848      CoinMemcpyN(columnActivityWork_,numberColumns_,savedSolution_);
4849      // save extra stuff
4850      int dummy;
4851      matrix_->generalExpanded(this,5,dummy);
4852    }
4853    if (weightsSaved) {
4854      // restore weights (if saved) - also recompute infeasibility list
4855      if (!reallyBadProblems&&(largestPrimalError_<100.0||numberPivots>10)) {
4856        if (tentativeStatus>-3) 
4857          dualRowPivot_->saveWeights(this,(type <2) ? 2 : 4);
4858        else
4859          dualRowPivot_->saveWeights(this,3);
4860      } else {
4861        // reset weights or scale back
4862        dualRowPivot_->saveWeights(this,6);
4863      }
4864    } else if (weightsSaved2&&numberPrimalInfeasibilities_) {
4865      dualRowPivot_->saveWeights(this,3);
4866    }
4867  }
4868  // see if cutoff reached
4869  double limit = 0.0;
4870  getDblParam(ClpDualObjectiveLimit, limit);
4871#if 0
4872  if(fabs(limit)<1.0e30&&objectiveValue()*optimizationDirection_>
4873     limit+100.0) {
4874    printf("lim %g obj %g %g - wo perturb %g sum dual %g\n",
4875           limit,objectiveValue_,objectiveValue(),computeInternalObjectiveValue(),sumDualInfeasibilities_);
4876  }
4877#endif
4878  if(fabs(limit)<1.0e30&&objectiveValue()*optimizationDirection_>
4879     limit&&!numberAtFakeBound()) {
4880    bool looksInfeasible = !numberDualInfeasibilities_;
4881    if (objectiveValue()*optimizationDirection_>limit+fabs(0.1*limit)+1.0e2*sumDualInfeasibilities_+1.0e4&&
4882        sumDualInfeasibilities_<largestDualError_&&numberIterations_>0.5*numberRows_+1000)
4883      looksInfeasible=true;
4884    if (looksInfeasible) {
4885      // Even if not perturbed internal costs may have changed
4886      // be careful
4887      if (true||numberIterations_) {
4888        if(computeInternalObjectiveValue()>limit) {
4889          problemStatus_=1;
4890          secondaryStatus_ = 1; // and say was on cutoff
4891        }
4892      } else {
4893        problemStatus_=1;
4894        secondaryStatus_ = 1; // and say was on cutoff
4895      }
4896    }
4897  }
4898  // If we are in trouble and in branch and bound give up
4899  if ((specialOptions_&1024)!=0) {
4900    int looksBad=0;
4901    if (largestPrimalError_*largestDualError_>1.0e2) {
4902      looksBad=1;
4903    } else if (largestPrimalError_>1.0e-2
4904        &&objectiveValue_>CoinMin(1.0e15,1.0e3*limit)) {
4905      looksBad=2;
4906    }
4907    if (looksBad) {
4908      if (factorization_->pivotTolerance()<0.9) {
4909        // up tolerance
4910        factorization_->pivotTolerance(CoinMin(factorization_->pivotTolerance()*1.05+0.02,0.91));
4911      } else if (numberIterations_>10000) {
4912        if (handler_->logLevel()>2)
4913          printf("bad dual - saying infeasible %d\n",looksBad);
4914        problemStatus_=1;
4915        secondaryStatus_ = 1; // and say was on cutoff
4916      } else if (largestPrimalError_>1.0e5) {
4917        {
4918          int iBigB=-1;
4919          double bigB=0.0;
4920          int iBigN=-1;
4921          double bigN=0.0;
4922          for (int i=0;i<numberRows_+numberColumns_;i++) {
4923            double value = fabs(solution_[i]);
4924            if (getStatus(i)==basic) {
4925              if (value>bigB) {
4926                bigB= value;
4927                iBigB=i;
4928              }
4929            } else {
4930              if (value>bigN) {
4931                bigN= value;
4932                iBigN=i;
4933              }
4934            }
4935          }
4936#ifdef CLP_INVESTIGATE
4937          if (bigB>1.0e8||bigN>1.0e8) {
4938            if (handler_->logLevel()>0)
4939              printf("it %d - basic %d %g, nonbasic %d %g\n",
4940                     numberIterations_,iBigB,bigB,iBigN,bigN);
4941          }
4942#endif
4943        }
4944#if COIN_DEVELOP!=2
4945        if (handler_->logLevel()>2)
4946#endif
4947          printf("bad dual - going to primal %d %g\n",looksBad,largestPrimalError_);
4948        allSlackBasis(true);
4949        problemStatus_=10;
4950      }
4951    }
4952  }
4953  if (problemStatus_<0&&!changeMade_) {
4954    problemStatus_=4; // unknown
4955  }
4956  lastGoodIteration_ = numberIterations_;
4957  if (numberIterations_>lastBadIteration_+100)
4958    moreSpecialOptions_ &= ~16; // clear check accuracy flag
4959  if (problemStatus_<0) {
4960    sumDualInfeasibilities_=realDualInfeasibilities; // back to say be careful
4961    if (sumDualInfeasibilities_)
4962      numberDualInfeasibilities_=1;
4963  }
4964#ifdef CLP_REPORT_PROGRESS
4965  if (ixxxxxx>ixxyyyy-3) {
4966    printf("objectiveValue_ %g\n",objectiveValue_);
4967    handler_->setLogLevel(63);
4968    int nTotal = numberColumns_+numberRows_;
4969    double newObj=0.0;
4970    for (int i=0;i<nTotal;i++) {
4971      if (solution_[i])
4972        newObj += solution_[i]*cost_[i];
4973    }
4974    printf("xxx obj %g\n",newObj);
4975    // for now - recompute all
4976    gutsOfSolution(NULL,NULL);
4977    newObj=0.0;
4978    for (int i=0;i<nTotal;i++) {
4979      if (solution_[i])
4980        newObj += solution_[i]*cost_[i];
4981    }
4982    printf("yyy obj %g %g\n",newObj,objectiveValue_);
4983    progress_.modifyObjective(objectiveValue_
4984                              -bestPossibleImprovement_);
4985  }
4986#endif
4987#if 1
4988  double thisObj = progress_.lastObjective(0);
4989  double lastObj = progress_.lastObjective(1);
4990  if (lastObj>thisObj+1.0e-4*CoinMax(fabs(thisObj),fabs(lastObj))+1.0e-4
4991      &&givenDuals==NULL&&firstFree_<0) {
4992    int maxFactor = factorization_->maximumPivots();
4993    if (maxFactor>10) {
4994      if (forceFactorization_<0)
4995        forceFactorization_= maxFactor;
4996      forceFactorization_ = CoinMax(1,(forceFactorization_>>1));
4997      //printf("Reducing factorization frequency\n");
4998    } 
4999  }
5000#endif
5001  // Allow matrices to be sorted etc
5002  int fake=-999; // signal sort
5003  matrix_->correctSequence(this,fake,fake);
5004  if (alphaAccuracy_>0.0)
5005      alphaAccuracy_=1.0;
5006}
5007/* While updateDualsInDual sees what effect is of flip
5008   this does actual flipping.
5009   If change >0.0 then value in array >0.0 => from lower to upper
5010*/
5011void 
5012ClpSimplexDual::flipBounds(CoinIndexedVector * rowArray,
5013                  CoinIndexedVector * columnArray,
5014                  double change)
5015{
5016  int number;
5017  int * which;
5018 
5019  int iSection;
5020
5021  for (iSection=0;iSection<2;iSection++) {
5022    int i;
5023    double * solution = solutionRegion(iSection);
5024    double * lower = lowerRegion(iSection);
5025    double * upper = upperRegion(iSection);
5026    int addSequence;
5027    if (!iSection) {
5028      number = rowArray->getNumElements();
5029      which = rowArray->getIndices();
5030      addSequence = numberColumns_;
5031    } else {
5032      number = columnArray->getNumElements();
5033      which = columnArray->getIndices();
5034      addSequence = 0;
5035    }
5036   
5037    for (i=0;i<number;i++) {
5038      int iSequence = which[i];
5039      Status status = getStatus(iSequence+addSequence);
5040
5041      switch(status) {
5042
5043      case basic:
5044      case isFree:
5045      case superBasic:
5046      case ClpSimplex::isFixed:
5047        break;
5048      case atUpperBound:
5049        // to lower bound
5050        setStatus(iSequence+addSequence,atLowerBound);
5051        solution[iSequence] = lower[iSequence];
5052        break;
5053      case atLowerBound:
5054        // to upper bound
5055        setStatus(iSequence+addSequence,atUpperBound);
5056        solution[iSequence] = upper[iSequence];
5057        break;
5058      }
5059    }
5060  }
5061  rowArray->setNumElements(0);
5062  columnArray->setNumElements(0);
5063}
5064// Restores bound to original bound
5065void 
5066ClpSimplexDual::originalBound( int iSequence)
5067{
5068  if (getFakeBound(iSequence)!=noFake) {
5069    numberFake_--;;
5070    setFakeBound(iSequence,noFake);
5071    if (iSequence>=numberColumns_) {
5072      // rows
5073      int iRow = iSequence-numberColumns_;
5074      rowLowerWork_[iRow]=rowLower_[iRow];
5075      rowUpperWork_[iRow]=rowUpper_[iRow];
5076      if (rowScale_) {
5077        if (rowLowerWork_[iRow]>-1.0e50)
5078          rowLowerWork_[iRow] *= rowScale_[iRow]*rhsScale_;
5079        if (rowUpperWork_[iRow]<1.0e50)
5080          rowUpperWork_[iRow] *= rowScale_[iRow]*rhsScale_;
5081      } else if (rhsScale_!=1.0) {
5082        if (rowLowerWork_[iRow]>-1.0e50)
5083          rowLowerWork_[iRow] *= rhsScale_;
5084        if (rowUpperWork_[iRow]<1.0e50)
5085          rowUpperWork_[iRow] *= rhsScale_;
5086      }
5087    } else {
5088      // columns
5089      columnLowerWork_[iSequence]=columnLower_[iSequence];
5090      columnUpperWork_[iSequence]=columnUpper_[iSequence];
5091      if (rowScale_) {
5092        double multiplier = 1.0*inverseColumnScale_[iSequence];
5093        if (columnLowerWork_[iSequence]>-1.0e50)
5094          columnLowerWork_[iSequence] *= multiplier*rhsScale_;
5095        if (columnUpperWork_[iSequence]<1.0e50)
5096          columnUpperWork_[iSequence] *= multiplier*rhsScale_;
5097      } else if (rhsScale_!=1.0) {
5098        if (columnLowerWork_[iSequence]>-1.0e50)
5099          columnLowerWork_[iSequence] *= rhsScale_;
5100        if (columnUpperWork_[iSequence]<1.0e50)
5101          columnUpperWork_[iSequence] *= rhsScale_;
5102      }
5103    }
5104  }
5105}
5106/* As changeBounds but just changes new bounds for a single variable.
5107   Returns true if change */
5108bool 
5109ClpSimplexDual::changeBound( int iSequence)
5110{
5111  // old values
5112  double oldLower=lower_[iSequence];
5113  double oldUpper=upper_[iSequence];
5114  double value=solution_[iSequence];
5115  bool modified=false;
5116  originalBound(iSequence);
5117  // original values
5118  double lowerValue=lower_[iSequence];
5119  double upperValue=upper_[iSequence];
5120  // back to altered values
5121  lower_[iSequence] = oldLower;
5122  upper_[iSequence] = oldUpper;
5123  assert (getFakeBound(iSequence)==noFake);
5124  //if (getFakeBound(iSequence)!=noFake)
5125  //numberFake_--;;
5126  if (value==oldLower) {
5127    if (upperValue > oldLower + dualBound_) {
5128      upper_[iSequence]=oldLower+dualBound_;
5129      setFakeBound(iSequence,upperFake);
5130      modified=true;
5131      numberFake_++;
5132    }
5133  } else if (value==oldUpper) {
5134    if (lowerValue < oldUpper - dualBound_) {
5135      lower_[iSequence]=oldUpper-dualBound_;
5136      setFakeBound(iSequence,lowerFake);
5137      modified=true;
5138      numberFake_++;
5139    }
5140  } else {
5141    assert(value==oldLower||value==oldUpper);
5142  }
5143  return modified;
5144}
5145// Perturbs problem
5146int 
5147ClpSimplexDual::perturb()
5148{
5149  if (perturbation_>100)
5150    return 0; //perturbed already
5151  if (perturbation_==100)
5152    perturbation_=50; // treat as normal
5153  int savePerturbation = perturbation_;
5154  bool modifyRowCosts=false;
5155  // dual perturbation
5156  double perturbation=1.0e-20;
5157  // maximum fraction of cost to perturb
5158  double maximumFraction = 1.0e-5;
5159  double constantPerturbation = 100.0*dualTolerance_;
5160  int maxLength=0;
5161  int minLength=numberRows_;
5162  double averageCost = 0.0;
5163#if 0
5164  // look at element range
5165  double smallestNegative;
5166  double largestNegative;
5167  double smallestPositive;
5168  double largestPositive;
5169  matrix_->rangeOfElements(smallestNegative, largestNegative,
5170                           smallestPositive, largestPositive);
5171  smallestPositive = CoinMin(fabs(smallestNegative),smallestPositive);
5172  largestPositive = CoinMax(fabs(largestNegative),largestPositive);
5173  double elementRatio = largestPositive/smallestPositive;
5174#endif
5175  int numberNonZero=0;
5176  if (!numberIterations_&&perturbation_>=50) {
5177    // See if we need to perturb
5178    double * sort = new double[numberColumns_];
5179    // Use objective BEFORE scaling
5180    const double * obj = objective();
5181    int i;
5182    for (i=0;i<numberColumns_;i++) {
5183      double value = fabs(obj[i]);
5184      sort[i]=value;
5185      averageCost += value;
5186      if (value)
5187        numberNonZero++;
5188    }
5189    if (numberNonZero)
5190      averageCost /= static_cast<double> (numberNonZero);
5191    else
5192      averageCost = 1.0;
5193    std::sort(sort,sort+numberColumns_);
5194    int number=1;
5195    double last = sort[0];
5196    for (i=1;i<numberColumns_;i++) {
5197      if (last!=sort[i])
5198        number++;
5199      last=sort[i];
5200    }
5201    delete [] sort;
5202    if (!numberNonZero)
5203      return 1; // safer to use primal
5204#if 0
5205    printf("nnz %d percent %d",number,(number*100)/numberColumns_);
5206    if (number*4>numberColumns_)
5207      printf(" - Would not perturb\n");
5208    else
5209      printf(" - Would perturb\n");
5210    //exit(0);
5211#endif
5212    //printf("ratio number diff costs %g, element ratio %g\n",((double)number)/((double) numberColumns_),
5213    //                                                                elementRatio);
5214    //number=0;
5215    //if (number*4>numberColumns_||elementRatio>1.0e12) {
5216    if (number*4>numberColumns_) {
5217      perturbation_=100;
5218      return 0; // good enough
5219    }
5220  }
5221  int iColumn;
5222  const int * columnLength = matrix_->getVectorLengths();
5223  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
5224    if (columnLowerWork_[iColumn]<columnUpperWork_[iColumn]) {
5225      int length = columnLength[iColumn];
5226      if (length>2) {
5227        maxLength = CoinMax(maxLength,length);
5228        minLength = CoinMin(minLength,length);
5229      }
5230    }
5231  }
5232  // If > 70 then do rows
5233  if (perturbation_>=70) {
5234    modifyRowCosts=true;
5235    perturbation_ -= 20;
5236    printf("Row costs modified, ");
5237  }
5238  bool uniformChange=false;
5239  bool inCbcOrOther = (specialOptions_&0x03000000)!=0;
5240  if (perturbation_>50) {
5241    // Experiment
5242    // maximumFraction could be 1.0e-10 to 1.0
5243    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};
5244    int whichOne = perturbation_-51;
5245    //if (inCbcOrOther&&whichOne>0)
5246    //whichOne--;
5247    maximumFraction = m[CoinMin(whichOne,10)];
5248  } else if (inCbcOrOther) {
5249    //maximumFraction = 1.0e-6;
5250  }
5251  int iRow;
5252  double smallestNonZero=1.0e100;
5253  numberNonZero=0;
5254  if (perturbation_>=50) {
5255    perturbation = 1.0e-8;
5256    bool allSame=true;
5257    double lastValue=0.0;
5258    for (iRow=0;iRow<numberRows_;iRow++) {
5259      double lo = rowLowerWork_[iRow];
5260      double up = rowUpperWork_[iRow];
5261      if (lo<up) {
5262        double value = fabs(rowObjectiveWork_[iRow]);
5263        perturbation = CoinMax(perturbation,value);
5264        if (value) {
5265          modifyRowCosts=true;
5266          smallestNonZero = CoinMin(smallestNonZero,value);
5267        }
5268      } 
5269      if (lo&&lo>-1.0e10) {
5270        numberNonZero++;
5271        lo=fabs(lo);
5272        if (!lastValue) 
5273          lastValue=lo;
5274        else if (fabs(lo-lastValue)>1.0e-7)
5275          allSame=false;
5276      }
5277      if (up&&up<1.0e10) {
5278        numberNonZero++;
5279        up=fabs(up);
5280        if (!lastValue) 
5281          lastValue=up;
5282        else if (fabs(up-lastValue)>1.0e-7)
5283          allSame=false;
5284      }
5285    }
5286    double lastValue2=0.0;
5287    for (iColumn=0;iColumn<numberColumns_;iColumn++) { 
5288      double lo = columnLowerWork_[iColumn];
5289      double up = columnUpperWork_[iColumn];
5290      if (lo<up) {
5291        double value = 
5292          fabs(objectiveWork_[iColumn]);
5293        perturbation = CoinMax(perturbation,value);
5294        if (value) {
5295          smallestNonZero = CoinMin(smallestNonZero,value);
5296        }
5297      }
5298      if (lo&&lo>-1.0e10) {
5299        //numberNonZero++;
5300        lo=fabs(lo);
5301        if (!lastValue2) 
5302          lastValue2=lo;
5303        else if (fabs(lo-lastValue2)>1.0e-7)
5304          allSame=false;
5305      }
5306      if (up&&up<1.0e10) {
5307        //numberNonZero++;
5308        up=fabs(up);
5309        if (!lastValue2) 
5310          lastValue2=up;
5311        else if (fabs(up-lastValue2)>1.0e-7)
5312          allSame=false;
5313      }
5314    }
5315    if (allSame) {
5316      // Check elements
5317      double smallestNegative;
5318      double largestNegative;
5319      double smallestPositive;
5320      double largestPositive;
5321      matrix_->rangeOfElements(smallestNegative,largestNegative,
5322                               smallestPositive,largestPositive);
5323      if (smallestNegative==largestNegative&&
5324          smallestPositive==largestPositive) {
5325        // Really hit perturbation
5326        double adjust = CoinMin(100.0*maximumFraction,1.0e-3*CoinMax(lastValue,lastValue2));
5327        maximumFraction = CoinMax(adjust,maximumFraction);
5328      }
5329    }
5330    perturbation = CoinMin(perturbation,smallestNonZero/maximumFraction);
5331  } else {
5332    // user is in charge
5333    maximumFraction = 1.0e-1;
5334    // but some experiments
5335    if (perturbation_<=-900) {
5336      modifyRowCosts=true;
5337      perturbation_ += 1000;
5338      printf("Row costs modified, ");
5339    }
5340    if (perturbation_<=-10) {
5341      perturbation_ += 10; 
5342      maximumFraction = 1.0;
5343      if ((-perturbation_)%100>=10) {
5344        uniformChange=true;
5345        perturbation_+=20;
5346      }
5347      while (perturbation_<-10) {
5348        perturbation_ += 100;
5349        maximumFraction *= 1.0e-1;
5350      }
5351    }
5352    perturbation = pow(10.0,perturbation_);
5353  }
5354  double largestZero=0.0;
5355  double largest=0.0;
5356  double largestPerCent=0.0;
5357  // modify costs
5358  bool printOut=(handler_->logLevel()==63);
5359  printOut=false;
5360  //assert (!modifyRowCosts);
5361  modifyRowCosts=false;
5362  if (modifyRowCosts) {
5363    for (iRow=0;iRow<numberRows_;iRow++) {
5364      if (rowLowerWork_[iRow]<rowUpperWork_[iRow]) {
5365        double value = perturbation;
5366        double currentValue = rowObjectiveWork_[iRow];
5367        value = CoinMin(value,maximumFraction*(fabs(currentValue)+1.0e-1*perturbation+1.0e-3));
5368        if (rowLowerWork_[iRow]>-largeValue_) {
5369          if (fabs(rowLowerWork_[iRow])<fabs(rowUpperWork_[iRow])) 
5370            value *= randomNumberGenerator_.randomDouble();
5371          else
5372            value *= -randomNumberGenerator_.randomDouble();
5373        } else if (rowUpperWork_[iRow]<largeValue_) {
5374          value *= -randomNumberGenerator_.randomDouble();
5375        } else {
5376          value=0.0;
5377        }
5378        if (currentValue) {
5379          largest = CoinMax(largest,fabs(value));
5380          if (fabs(value)>fabs(currentValue)*largestPerCent)
5381            largestPerCent=fabs(value/currentValue);
5382        } else {
5383          largestZero=CoinMax(largestZero,fabs(value));
5384        }
5385        if (printOut)
5386          printf("row %d cost %g change %g\n",iRow,rowObjectiveWork_[iRow],value);
5387        rowObjectiveWork_[iRow] += value;
5388      }
5389    }
5390  }
5391  // 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};
5392  // 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};
5393  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,20.0,50.0,100.0,120.0,130.0,140.0,200.0};
5395  double extraWeight=10.0;
5396  // Scale back if wanted
5397  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};
5398  if (constantPerturbation<99.0*dualTolerance_) {
5399    perturbation *= 0.1;
5400    extraWeight=0.5;
5401    memcpy(weight,weight2,sizeof(weight2));
5402  }
5403  // adjust weights if all columns long
5404  double factor=1.0;
5405  if (maxLength) {
5406    factor = 3.0/static_cast<double> (minLength);
5407  }
5408  // Make variables with more elements more expensive
5409  const double m1 = 0.5;
5410  double smallestAllowed = CoinMin(1.0e-2*dualTolerance_,maximumFraction);
5411  double largestAllowed = CoinMax(1.0e3*dualTolerance_,maximumFraction*averageCost);
5412  // smaller if in BAB
5413  //if (inCbcOrOther)
5414  //largestAllowed=CoinMin(largestAllowed,1.0e-5);
5415  //smallestAllowed = CoinMin(smallestAllowed,0.1*largestAllowed);
5416#define SAVE_PERT
5417#ifdef SAVE_PERT
5418  if (2*numberColumns_>maximumPerturbationSize_) {
5419    delete [] perturbationArray_;
5420    maximumPerturbationSize_ = 2* numberColumns_;
5421    perturbationArray_ = new double [maximumPerturbationSize_];
5422    for (iColumn=0;iColumn<maximumPerturbationSize_;iColumn++) {
5423      perturbationArray_[iColumn] = randomNumberGenerator_.randomDouble();
5424    }
5425  }
5426#endif
5427  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
5428    if (columnLowerWork_[iColumn]<columnUpperWork_[iColumn]&&getStatus(iColumn)!=basic) {
5429      double value = perturbation;
5430      double currentValue = objectiveWork_[iColumn];
5431      value = CoinMin(value,constantPerturbation+maximumFraction*(fabs(currentValue)+1.0e-1*perturbation+1.0e-8));
5432      //value = CoinMin(value,constantPerturbation;+maximumFraction*fabs(currentValue));
5433      double value2 = constantPerturbation+1.0e-1*smallestNonZero;
5434      if (uniformChange) {
5435        value = maximumFraction;
5436        value2=maximumFraction;
5437      }
5438      if (columnLowerWork_[iColumn]>-largeValue_) {
5439        if (fabs(columnLowerWork_[iColumn])<
5440            fabs(columnUpperWork_[iColumn])) {
5441#ifndef SAVE_PERT
5442          value *= (1.0-m1 +m1*randomNumberGenerator_.randomDouble());
5443          value2 *= (1.0-m1 +m1*randomNumberGenerator_.randomDouble());
5444#else
5445          value *= (1.0-m1 +m1*perturbationArray_[2*iColumn]);
5446          value2 *= (1.0-m1 +m1*perturbationArray_[2*iColumn+1]);
5447#endif
5448        } else {
5449          //value *= -(1.0-m1+m1*randomNumberGenerator_.randomDouble());
5450          //value2 *= -(1.0-m1+m1*randomNumberGenerator_.randomDouble());
5451          value=0.0;
5452        }
5453      } else if (columnUpperWork_[iColumn]<largeValue_) {
5454#ifndef SAVE_PERT
5455        value *= -(1.0-m1+m1*randomNumberGenerator_.randomDouble());
5456        value2 *= -(1.0-m1+m1*randomNumberGenerator_.randomDouble());
5457#else
5458        value *= -(1.0-m1+m1*perturbationArray_[2*iColumn]);
5459        value2 *= -(1.0-m1+m1*perturbationArray_[2*iColumn+1]);
5460#endif
5461      } else {
5462        value=0.0;
5463      }
5464      if (value) {
5465        int length = columnLength[iColumn];
5466        if (length>3) {
5467          length = static_cast<int> (static_cast<double> (length) * factor);
5468          length = CoinMax(3,length);
5469        }
5470        double multiplier;
5471#if 1
5472        if (length<10)
5473          multiplier=weight[length];
5474        else
5475          multiplier = weight[10];
5476#else
5477        if (length<10)
5478          multiplier=weight[length];
5479        else
5480          multiplier = weight[10]+extraWeight*(length-10);
5481        multiplier *= 0.5;
5482#endif
5483        value *= multiplier;
5484        value = CoinMin(value,value2);
5485        if (savePerturbation<50||savePerturbation>60) {
5486          if (fabs(value)<=dualTolerance_)
5487            value=0.0;
5488        } else if (value) {
5489          // get in range
5490          if (fabs(value)<=smallestAllowed) {
5491            value *= 10.0;
5492            while (fabs(value)<=smallestAllowed) 
5493              value *= 10.0;
5494          } else if (fabs(value)>largestAllowed) {
5495            value *= 0.1;
5496            while (fabs(value)>largestAllowed) 
5497              value *= 0.1;
5498          }
5499        }
5500        if (currentValue) {
5501          largest = CoinMax(largest,fabs(value));
5502          if (fabs(value)>fabs(currentValue)*largestPerCent)
5503            largestPerCent=fabs(value/currentValue);
5504        } else {
5505          largestZero=CoinMax(largestZero,fabs(value));
5506        }
5507        // but negative if at ub
5508        if (getStatus(iColumn)==atUpperBound)
5509          value = -value;
5510        if (printOut)
5511          printf("col %d cost %g change %g\n",iColumn,objectiveWork_[iColumn],value);
5512        objectiveWork_[iColumn] += value;
5513      }
5514    }
5515  }
5516  handler_->message(CLP_SIMPLEX_PERTURB,messages_)
5517    <<100.0*maximumFraction<<perturbation<<largest<<100.0*largestPerCent<<largestZero
5518    <<CoinMessageEol;
5519  // and zero changes
5520  //int nTotal = numberRows_+numberColumns_;
5521  //CoinZeroN(cost_+nTotal,nTotal);
5522  // say perturbed
5523  perturbation_=101;
5524  return 0;
5525}
5526/* For strong branching.  On input lower and upper are new bounds
5527   while on output they are change in objective function values
5528   (>1.0e50 infeasible).
5529   Return code is 0 if nothing interesting, -1 if infeasible both
5530   ways and +1 if infeasible one way (check values to see which one(s))
5531   Returns -2 if bad factorization
5532*/
5533int ClpSimplexDual::strongBranching(int numberVariables,const int * variables,
5534                                    double * newLower, double * newUpper,
5535                                    double ** outputSolution,
5536                                    int * outputStatus, int * outputIterations,
5537                                    bool stopOnFirstInfeasible,
5538                                    bool alwaysFinish,
5539                                    int startFinishOptions)
5540{
5541  int i;
5542  int returnCode=0;
5543  double saveObjectiveValue = objectiveValue_;
5544  algorithm_ = -1;
5545
5546  //scaling(false);
5547
5548  // put in standard form (and make row copy)
5549  // create modifiable copies of model rim and do optional scaling
5550  createRim(7+8+16+32,true,startFinishOptions);
5551
5552  // change newLower and newUpper if scaled
5553
5554  // Do initial factorization
5555  // and set certain stuff
5556  // We can either set increasing rows so ...IsBasic gives pivot row
5557  // or we can just increment iBasic one by one
5558  // for now let ...iBasic give pivot row
5559  int useFactorization=false;
5560  if ((startFinishOptions&2)!=0&&(whatsChanged_&(2+512))==2+512)
5561    useFactorization=true; // Keep factorization if possible
5562  // switch off factorization if bad
5563  if (pivotVariable_[0]<0)
5564    useFactorization=false;
5565  if (!useFactorization||factorization_->numberRows()!=numberRows_) {
5566    useFactorization = false;
5567    factorization_->setDefaultValues();
5568
5569    int factorizationStatus = internalFactorize(0);
5570    if (factorizationStatus<0) {
5571      // some error
5572      // we should either debug or ignore
5573#ifndef NDEBUG
5574      printf("***** ClpDual strong branching factorization error - debug\n");
5575#endif
5576      return -2;
5577    } else if (factorizationStatus&&factorizationStatus<=numberRows_) {
5578      handler_->message(CLP_SINGULARITIES,messages_)
5579        <<factorizationStatus
5580        <<CoinMessageEol;
5581    }
5582  }
5583  // save stuff
5584  ClpFactorization saveFactorization(*factorization_);
5585  // save basis and solution
5586  double * saveSolution = new double[numberRows_+numberColumns_];
5587  CoinMemcpyN(solution_,
5588         numberRows_+numberColumns_,saveSolution);
5589  unsigned char * saveStatus = 
5590    new unsigned char [numberRows_+numberColumns_];
5591  CoinMemcpyN(status_,numberColumns_+numberRows_,saveStatus);
5592  // save bounds as createRim makes clean copies
5593  double * saveLower = new double[numberRows_+numberColumns_];
5594  CoinMemcpyN(lower_,
5595         numberRows_+numberColumns_,saveLower);
5596  double * saveUpper = new double[numberRows_+numberColumns_];
5597  CoinMemcpyN(upper_,
5598         numberRows_+numberColumns_,saveUpper);
5599  double * saveObjective = new double[numberRows_+numberColumns_];
5600  CoinMemcpyN(cost_,
5601         numberRows_+numberColumns_,saveObjective);
5602  int * savePivot = new int [numberRows_];
5603  CoinMemcpyN(pivotVariable_, numberRows_,savePivot);
5604  // need to save/restore weights.
5605
5606  int iSolution = 0;
5607  for (i=0;i<numberVariables;i++) {
5608    int iColumn = variables[i];
5609    double objectiveChange;
5610    double saveBound;
5611   
5612    // try down
5613   
5614    saveBound = columnUpper_[iColumn];
5615    // external view - in case really getting optimal
5616    columnUpper_[iColumn] = newUpper[i];
5617    assert (inverseColumnScale_);
5618    if (scalingFlag_<=0) 
5619      upper_[iColumn] = newUpper[i]*rhsScale_;
5620    else 
5621      upper_[iColumn] = (newUpper[i]*inverseColumnScale_[iColumn])*rhsScale_; // scale
5622    // Get fake bounds correctly
5623    double changeCost;
5624    changeBounds(3,NULL,changeCost);
5625    // Start of fast iterations
5626    int status = fastDual(alwaysFinish);
5627    CoinAssert (problemStatus_||objectiveValue_<1.0e50);
5628    // make sure plausible
5629    double obj = CoinMax(objectiveValue_,saveObjectiveValue);
5630    if (status&&problemStatus_!=3) {
5631      // not finished - might be optimal
5632      checkPrimalSolution(rowActivityWork_,columnActivityWork_);
5633      double limit = 0.0;
5634      getDblParam(ClpDualObjectiveLimit, limit);
5635      if (!numberPrimalInfeasibilities_&&obj<limit) { 
5636        problemStatus_=0;
5637      } 
5638      status=problemStatus_;
5639    }
5640    if (status||(problemStatus_==0&&!isDualObjectiveLimitReached())) {
5641      objectiveChange = obj-saveObjectiveValue;
5642    } else {
5643      objectiveChange = 1.0e100;
5644      status=1;
5645    }
5646    if (problemStatus_==3)
5647      status=2;
5648
5649    if (scalingFlag_<=0) {
5650      CoinMemcpyN(solution_,numberColumns_,outputSolution[iSolution]);
5651    } else {
5652      int j;
5653      double * sol = outputSolution[iSolution];
5654      for (j=0;j<numberColumns_;j++) 
5655        sol[j] = solution_[j]*columnScale_[j];
5656    }
5657    outputStatus[iSolution]=status;
5658    outputIterations[iSolution]=numberIterations_;
5659    iSolution++;
5660    // restore
5661    CoinMemcpyN(saveSolution,
5662            numberRows_+numberColumns_,solution_);
5663    CoinMemcpyN(saveStatus,numberColumns_+numberRows_,status_);
5664    CoinMemcpyN(saveLower,
5665            numberRows_+numberColumns_,lower_);
5666    CoinMemcpyN(saveUpper,
5667            numberRows_+numberColumns_,upper_);
5668    CoinMemcpyN(saveObjective,
5669            numberRows_+numberColumns_,cost_);
5670    columnUpper_[iColumn] = saveBound;
5671    CoinMemcpyN(savePivot, numberRows_,pivotVariable_);
5672    delete factorization_;
5673    factorization_ = new ClpFactorization(saveFactorization,numberRows_);
5674
5675    newUpper[i]=objectiveChange;
5676#ifdef CLP_DEBUG
5677    printf("down on %d costs %g\n",iColumn,objectiveChange);
5678#endif
5679         
5680    // try up
5681   
5682    saveBound = columnLower_[iColumn];
5683    // external view - in case really getting optimal
5684    columnLower_[iColumn] = newLower[i];
5685    assert (inverseColumnScale_);
5686    if (scalingFlag_<=0) 
5687      lower_[iColumn] = newLower[i]*rhsScale_;
5688    else 
5689      lower_[iColumn] = (newLower[i]*inverseColumnScale_[iColumn])*rhsScale_; // scale
5690    // Get fake bounds correctly
5691    resetFakeBounds(1);
5692    // Start of fast iterations
5693    status = fastDual(alwaysFinish);
5694    // make sure plausible
5695    obj = CoinMax(objectiveValue_,saveObjectiveValue);
5696    if (status&&problemStatus_!=3) {
5697      // not finished - might be optimal
5698      checkPrimalSolution(rowActivityWork_,columnActivityWork_);
5699      double limit = 0.0;
5700      getDblParam(ClpDualObjectiveLimit, limit);
5701      if (!numberPrimalInfeasibilities_&&obj< limit) { 
5702        problemStatus_=0;
5703      } 
5704      status=problemStatus_;
5705    }
5706    if (status||(problemStatus_==0&&!isDualObjectiveLimitReached())) {
5707      objectiveChange = obj-saveObjectiveValue;
5708    } else {
5709      objectiveChange = 1.0e100;
5710      status=1;
5711    }
5712    if (problemStatus_==3)
5713      status=2;
5714    if (scalingFlag_<=0) {
5715      CoinMemcpyN(solution_,numberColumns_,outputSolution[iSolution]);
5716    } else {
5717      int j;
5718      double * sol = outputSolution[iSolution];
5719      for (j=0;j<numberColumns_;j++) 
5720        sol[j] = solution_[j]*columnScale_[j];
5721    }
5722    outputStatus[iSolution]=status;
5723    outputIterations[iSolution]=numberIterations_;
5724    iSolution++;
5725
5726    // restore
5727    CoinMemcpyN(saveSolution,
5728            numberRows_+numberColumns_,solution_);
5729    CoinMemcpyN(saveStatus,numberColumns_+numberRows_,status_);
5730    CoinMemcpyN(saveLower,
5731            numberRows_+numberColumns_,lower_);
5732    CoinMemcpyN(saveUpper,
5733            numberRows_+numberColumns_,upper_);
5734    CoinMemcpyN(saveObjective,
5735            numberRows_+numberColumns_,cost_);
5736    columnLower_[iColumn] = saveBound;
5737    CoinMemcpyN(savePivot, numberRows_,pivotVariable_);
5738    delete factorization_;
5739    factorization_ = new ClpFactorization(saveFactorization,numberRows_);
5740
5741    newLower[i]=objectiveChange;
5742#ifdef CLP_DEBUG
5743    printf("up on %d costs %g\n",iColumn,objectiveChange);
5744#endif
5745         
5746    /* Possibilities are:
5747       Both sides feasible - store
5748       Neither side feasible - set objective high and exit if desired
5749       One side feasible - change bounds and resolve
5750    */
5751    if (newUpper[i]<1.0e100) {
5752      if(newLower[i]<1.0e100) {
5753        // feasible - no action
5754      } else {
5755        // up feasible, down infeasible
5756        returnCode=1;
5757        if (stopOnFirstInfeasible)
5758          break;
5759      }
5760    } else {
5761      if(newLower[i]<1.0e100) {
5762        // down feasible, up infeasible
5763        returnCode=1;
5764        if (stopOnFirstInfeasible)
5765          break;
5766      } else {
5767        // neither side feasible
5768        returnCode=-1;
5769        break;
5770      }
5771    }
5772  }
5773  delete [] saveSolution;
5774  delete [] saveLower;
5775  delete [] saveUpper;
5776  delete [] saveObjective;
5777  delete [] saveStatus;
5778  delete [] savePivot;
5779  if ((startFinishOptions&1)==0) {
5780    deleteRim(1);
5781    whatsChanged_ &= ~0xffff;
5782  } else {
5783    // Original factorization will have been put back by last loop
5784    //delete factorization_;
5785    //factorization_ = new ClpFactorization(saveFactorization);
5786    deleteRim(0);
5787    // mark all as current
5788    whatsChanged_ = 0x3ffffff;
5789  }
5790  objectiveValue_ = saveObjectiveValue;
5791  return returnCode;
5792}
5793// treat no pivot as finished (unless interesting)
5794int