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

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

changes to make dual a nit more reliable

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