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

Last change on this file since 1321 was 1321, checked in by forrest, 13 years ago

out compiler warnings and stability improvements

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