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

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

bit more safety on fake bound handling

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 210.3 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4
5/* Notes on implementation of dual simplex algorithm.
6
7   When dual feasible:
8
9   If primal feasible, we are optimal.  Otherwise choose an infeasible
10   basic variable to leave basis (normally going to nearest bound) (B).  We
11   now need to find an incoming variable which will leave problem
12   dual feasible so we get the row of the tableau corresponding to
13   the basic variable (with the correct sign depending if basic variable
14   above or below feasibility region - as that affects whether reduced
15   cost on outgoing variable has to be positive or negative).
16
17   We now perform a ratio test to determine which incoming variable will
18   preserve dual feasibility (C).  If no variable found then problem
19   is infeasible (in primal sense).  If there is a variable, we then
20   perform pivot and repeat.  Trivial?
21
22   -------------------------------------------
23
24   A) How do we get dual feasible?  If all variables have bounds then
25   it is trivial to get feasible by putting non-basic variables to
26   correct bounds.  OSL did not have a phase 1/phase 2 approach but
27   instead effectively put fake bounds on variables and this is the
28   approach here, although I had hoped to make it cleaner.
29
30   If there is a weight of X on getting dual feasible:
31     Non-basic variables with negative reduced costs are put to
32     lesser of their upper bound and their lower bound + X.
33     Similarly, mutatis mutandis, for positive reduced costs.
34
35   Free variables should normally be in basis, otherwise I have
36   coding which may be able to come out (and may not be correct).
37
38   In OSL, this weight was changed heuristically, here at present
39   it is only increased if problem looks finished.  If problem is
40   feasible I check for unboundedness.  If not unbounded we
41   could play with going into primal.  As long as weights increase
42   any algorithm would be finite.
43   
44   B) Which outgoing variable to choose is a virtual base class.
45   For difficult problems steepest edge is preferred while for
46   very easy (large) problems we will need partial scan.
47
48   C) Sounds easy, but this is hardest part of algorithm.
49      1) Instead of stopping at first choice, we may be able
50      to flip that variable to other bound and if objective
51      still improving choose again.  These mini iterations can
52      increase speed by orders of magnitude but we may need to
53      go to more of a bucket choice of variable rather than looking
54      at them one by one (for speed).
55      2) Accuracy.  Reduced costs may be of wrong sign but less than
56      tolerance.  Pivoting on these makes objective go backwards.
57      OSL modified cost so a zero move was made, Gill et al
58      (in primal analogue) modified so a strictly positive move was
59      made.  It is not quite as neat in dual but that is what we
60      try and do.  The two problems are that re-factorizations can
61      change reduced costs above and below tolerances and that when
62      finished we need to reset costs and try again.
63      3) Degeneracy.  Gill et al helps but may not be enough.  We
64      may need more.  Also it can improve speed a lot if we perturb
65      the costs significantly. 
66
67  References:
68     Forrest and Goldfarb, Steepest-edge simplex algorithms for
69       linear programming - Mathematical Programming 1992
70     Forrest and Tomlin, Implementing the simplex method for
71       the Optimization Subroutine Library - IBM Systems Journal 1992
72     Gill, Murray, Saunders, Wright A Practical Anti-Cycling
73       Procedure for Linear and Nonlinear Programming SOL report 1988
74
75
76  TODO:
77 
78  a) Better recovery procedures.  At present I never check on forward
79     progress.  There is checkpoint/restart with reducing
80     re-factorization frequency, but this is only on singular
81     factorizations.
82  b) Fast methods for large easy problems (and also the option for
83     the code to automatically choose which method).
84  c) We need to be able to stop in various ways for OSI - this
85     is fairly easy.
86
87 */
88#ifdef COIN_DEVELOP
89#undef COIN_DEVELOP
90#define COIN_DEVELOP 2
91#endif
92
93#include "CoinPragma.hpp"
94
95#include <math.h>
96
97#include "CoinHelperFunctions.hpp"
98#include "ClpHelperFunctions.hpp"
99#include "ClpSimplexDual.hpp"
100#include "ClpEventHandler.hpp"
101#include "ClpFactorization.hpp"
102#include "CoinPackedMatrix.hpp"
103#include "CoinIndexedVector.hpp"
104#include "CoinFloatEqual.hpp"
105#include "ClpDualRowDantzig.hpp"
106#include "ClpMessage.hpp"
107#include "ClpLinearObjective.hpp"
108#include <cfloat>
109#include <cassert>
110#include <string>
111#include <stdio.h>
112#include <iostream>
113//#define CLP_DEBUG 1
114// To force to follow another run put logfile name here and define
115//#define FORCE_FOLLOW
116#ifdef FORCE_FOLLOW
117static FILE * fpFollow=NULL;
118static char * forceFile="old.log";
119static int force_in=-1;
120static int force_out=-1;
121static int force_iteration=0;
122#endif
123//#define VUB
124#ifdef VUB
125extern int * vub;
126extern int * toVub;
127extern int * nextDescendent;
128#endif
129#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  int dontFactorizePivots = dontFactorizePivots_;
3863  if (type==3) {
3864    type=1;
3865    dontFactorizePivots=1;
3866  }
3867  if (alphaAccuracy_<0.0||!numberPivots||alphaAccuracy_>1.0e4||numberPivots>20) {
3868    if (problemStatus_>-3||numberPivots>dontFactorizePivots) {
3869      // factorize
3870      // later on we will need to recover from singularities
3871      // also we could skip if first time
3872      // save dual weights
3873      dualRowPivot_->saveWeights(this,1);
3874      weightsSaved=true;
3875      if (type) {
3876        // is factorization okay?
3877        if (internalFactorize(1)) {
3878          // no - restore previous basis
3879          unflagVariables = false;
3880          assert (type==1);
3881          changeMade_++; // say something changed
3882          // Keep any flagged variables
3883          int i;
3884          for (i=0;i<numberRows_+numberColumns_;i++) {
3885            if (flagged(i))
3886              saveStatus_[i] |= 64; //say flagged
3887          }
3888          CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
3889          CoinMemcpyN(savedSolution_+numberColumns_ ,
3890                      numberRows_,rowActivityWork_);
3891          CoinMemcpyN(savedSolution_ ,
3892                      numberColumns_,columnActivityWork_);
3893          // restore extra stuff
3894          int dummy;
3895          matrix_->generalExpanded(this,6,dummy);
3896          // get correct bounds on all variables
3897          resetFakeBounds(1);
3898          // need to reject something
3899          char x = isColumn(sequenceOut_) ? 'C' :'R';
3900          handler_->message(CLP_SIMPLEX_FLAG,messages_)
3901            <<x<<sequenceWithin(sequenceOut_)
3902            <<CoinMessageEol;
3903#ifdef COIN_DEVELOP
3904          printf("flag d\n");
3905#endif
3906          setFlagged(sequenceOut_);
3907          progress_.clearBadTimes();
3908         
3909          // Go to safe
3910          factorization_->pivotTolerance(0.99);
3911          forceFactorization_=1; // a bit drastic but ..
3912          type = 2;
3913          //assert (internalFactorize(1)==0);
3914          if (internalFactorize(1)) {
3915            CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
3916            CoinMemcpyN(savedSolution_+numberColumns_ ,
3917                        numberRows_,rowActivityWork_);
3918            CoinMemcpyN(savedSolution_ ,
3919                        numberColumns_,columnActivityWork_);
3920            // restore extra stuff
3921            int dummy;
3922            matrix_->generalExpanded(this,6,dummy);
3923            // debug
3924            int returnCode = internalFactorize(1);
3925            while (returnCode) {
3926              // ouch
3927              // switch off dense
3928              int saveDense = factorization_->denseThreshold();
3929              factorization_->setDenseThreshold(0);
3930              // Go to safe
3931              factorization_->pivotTolerance(0.99);
3932              // make sure will do safe factorization
3933              pivotVariable_[0]=-1;
3934              returnCode=internalFactorize(2);
3935              factorization_->setDenseThreshold(saveDense);
3936            }
3937            // get correct bounds on all variables
3938            resetFakeBounds(1);
3939          }
3940        }
3941      }
3942      if (problemStatus_!=-4||numberPivots>10)
3943        problemStatus_=-3;
3944    }
3945  } else {
3946    //printf("testing with accuracy of %g and status of %d\n",alphaAccuracy_,problemStatus_);
3947    //count_alpha++;
3948    //if ((count_alpha%5000)==0)
3949    //printf("count alpha %d\n",count_alpha);
3950  }
3951  // at this stage status is -3 or -4 if looks infeasible
3952  // get primal and dual solutions
3953#if 0
3954  {
3955    int numberTotal = numberRows_+numberColumns_;
3956    double * saveSol = CoinCopyOfArray(solution_,numberTotal);
3957    double * saveDj = CoinCopyOfArray(dj_,numberTotal);
3958    double tolerance = type ? 1.0e-4 : 1.0e-8;
3959    // always if values pass
3960    double saveObj=objectiveValue_;
3961    double sumPrimal = sumPrimalInfeasibilities_;
3962    int numberPrimal = numberPrimalInfeasibilities_;
3963    double sumDual = sumDualInfeasibilities_;
3964    int numberDual = numberDualInfeasibilities_;
3965    gutsOfSolution(givenDuals,NULL);
3966    int j;
3967    double largestPrimal = tolerance;
3968    int iPrimal=-1;
3969    for (j=0;j<numberTotal;j++) {
3970      double difference = solution_[j]-saveSol[j];
3971      if (fabs(difference)>largestPrimal) {
3972        iPrimal=j;
3973        largestPrimal=fabs(difference);
3974      }
3975    }
3976    double largestDual = tolerance;
3977    int iDual=-1;
3978    for (j=0;j<numberTotal;j++) {
3979      double difference = dj_[j]-saveDj[j];
3980      if (fabs(difference)>largestDual&&upper_[j]>lower_[j]) {
3981        iDual=j;
3982        largestDual=fabs(difference);
3983      }
3984    }
3985    if (!type) {
3986      if (fabs(saveObj-objectiveValue_)>1.0e-5||
3987          numberPrimal!=numberPrimalInfeasibilities_||numberPrimal!=1||
3988          fabs(sumPrimal-sumPrimalInfeasibilities_)>1.0e-5||iPrimal>=0||
3989          numberDual!=numberDualInfeasibilities_||numberDual!=0||
3990          fabs(sumDual-sumDualInfeasibilities_)>1.0e-5||iDual>=0)
3991        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",
3992               type,numberIterations_,numberPivots,
3993               numberPrimal,numberPrimalInfeasibilities_,sumPrimal,sumPrimalInfeasibilities_,
3994               largestPrimal,iPrimal,
3995               numberDual,numberDualInfeasibilities_,sumDual,sumDualInfeasibilities_,
3996               largestDual,iDual,
3997               saveObj,objectiveValue_);
3998    } else {
3999      if (fabs(saveObj-objectiveValue_)>1.0e-5||
4000          numberPrimalInfeasibilities_||iPrimal>=0||
4001          numberDualInfeasibilities_||iDual>=0)
4002        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",
4003               type,numberIterations_,numberPivots,
4004               numberPrimal,numberPrimalInfeasibilities_,sumPrimal,sumPrimalInfeasibilities_,
4005               largestPrimal,iPrimal,
4006               numberDual,numberDualInfeasibilities_,sumDual,sumDualInfeasibilities_,
4007               largestDual,iDual,
4008               saveObj,objectiveValue_);
4009    }
4010    delete [] saveSol;
4011    delete [] saveDj;
4012  }
4013#else
4014  if (type||ifValuesPass)
4015    gutsOfSolution(givenDuals,NULL);
4016#endif
4017  // If bad accuracy treat as singular
4018  if ((largestPrimalError_>1.0e15||largestDualError_>1.0e15)&&numberIterations_) {
4019    // restore previous basis
4020    unflagVariables = false;
4021    changeMade_++; // say something changed
4022    // Keep any flagged variables
4023    int i;
4024    for (i=0;i<numberRows_+numberColumns_;i++) {
4025      if (flagged(i))
4026        saveStatus_[i] |= 64; //say flagged
4027    }
4028    CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
4029    CoinMemcpyN(savedSolution_+numberColumns_ ,
4030           numberRows_,rowActivityWork_);
4031    CoinMemcpyN(savedSolution_ ,
4032           numberColumns_,columnActivityWork_);
4033    // restore extra stuff
4034    int dummy;
4035    matrix_->generalExpanded(this,6,dummy);
4036    // get correct bounds on all variables
4037    resetFakeBounds(1);
4038    // need to reject something
4039    char x = isColumn(sequenceOut_) ? 'C' :'R';
4040    handler_->message(CLP_SIMPLEX_FLAG,messages_)
4041      <<x<<sequenceWithin(sequenceOut_)
4042      <<CoinMessageEol;
4043#ifdef COIN_DEVELOP
4044    printf("flag e\n");
4045#endif
4046    setFlagged(sequenceOut_);
4047    progress_.clearBadTimes();
4048   
4049    // Go to safer
4050    double newTolerance = CoinMin(1.1*factorization_->pivotTolerance(),0.99);
4051    factorization_->pivotTolerance(newTolerance);
4052    forceFactorization_=1; // a bit drastic but ..
4053    if (alphaAccuracy_!=-1.0)
4054      alphaAccuracy_=-2.0;
4055    type = 2;
4056    //assert (internalFactorize(1)==0);
4057    if (internalFactorize(1)) {
4058      CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
4059      CoinMemcpyN(savedSolution_+numberColumns_ ,
4060              numberRows_,rowActivityWork_);
4061      CoinMemcpyN(savedSolution_ ,
4062           numberColumns_,columnActivityWork_);
4063      // restore extra stuff
4064      int dummy;
4065      matrix_->generalExpanded(this,6,dummy);
4066      // debug
4067      int returnCode = internalFactorize(1);
4068      while (returnCode) {
4069        // ouch
4070        // switch off dense
4071        int saveDense = factorization_->denseThreshold();
4072        factorization_->setDenseThreshold(0);
4073        // Go to safe
4074        factorization_->pivotTolerance(0.99);
4075        // make sure will do safe factorization
4076        pivotVariable_[0]=-1;
4077        returnCode=internalFactorize(2);
4078        factorization_->setDenseThreshold(saveDense);
4079      }
4080      // get correct bounds on all variables
4081      resetFakeBounds(1);
4082    }
4083    // get primal and dual solutions
4084    gutsOfSolution(givenDuals,NULL);
4085  } else if (goodAccuracy()) {
4086    // Can reduce tolerance
4087    double newTolerance = CoinMax(0.99*factorization_->pivotTolerance(),saveData.pivotTolerance_);
4088    factorization_->pivotTolerance(newTolerance);
4089  }
4090  bool reallyBadProblems=false;
4091  // Double check infeasibility if no action
4092  if (progress_.lastIterationNumber(0)==numberIterations_) {
4093    if (dualRowPivot_->looksOptimal()) {
4094      numberPrimalInfeasibilities_ = 0;
4095      sumPrimalInfeasibilities_ = 0.0;
4096    }
4097#if 1
4098  } else {
4099    double thisObj = objectiveValue_-bestPossibleImprovement_;
4100#ifdef CLP_INVESTIGATE
4101    if (bestPossibleImprovement_)
4102      printf("obj %g add in %g -> %g\n",objectiveValue_,bestPossibleImprovement_,
4103             thisObj);
4104#endif
4105    double lastObj = progress_.lastObjective(0);
4106#ifndef NDEBUG
4107#ifdef COIN_DEVELOP
4108    resetFakeBounds(-1);
4109#endif
4110#endif
4111#ifdef CLP_REPORT_PROGRESS
4112    ixxxxxx++;
4113    if (ixxxxxx>=ixxyyyy-4&&ixxxxxx<=ixxyyyy) {
4114      char temp[20];
4115      sprintf(temp,"sol%d.out",ixxxxxx);
4116      printf("sol%d.out\n",ixxxxxx);
4117      FILE * fp = fopen(temp,"w");
4118      int nTotal=numberRows_+numberColumns_;
4119      for (int i=0;i<nTotal;i++)
4120        fprintf(fp,"%d %d %g %g %g %g %g\n",
4121                i,status_[i],lower_[i],solution_[i],upper_[i],cost_[i],dj_[i]);
4122      fclose(fp);
4123    }
4124#endif
4125    if(!ifValuesPass&&firstFree_<0) {
4126      double testTol = 5.0e-3;
4127      if (progress_.timesFlagged()>10) {
4128        testTol *= pow(2.0,progress_.timesFlagged()-8);
4129      } else if (progress_.timesFlagged()>5) {
4130        testTol *=5.0;
4131      }
4132      if (lastObj>thisObj+
4133          testTol*(fabs(thisObj)+fabs(lastObj))+testTol) {
4134        int maxFactor = factorization_->maximumPivots();
4135        if ((specialOptions_&16384)==0) {
4136          if (progress_.timesFlagged()>10) 
4137            progress_.incrementReallyBadTimes();
4138          if (maxFactor>10-9) {
4139#ifdef COIN_DEVELOP
4140            printf("lastobj %g thisobj %g\n",lastObj,thisObj);
4141#endif
4142            //if (forceFactorization_<0)
4143            //forceFactorization_= maxFactor;
4144            //forceFactorization_ = CoinMax(1,(forceFactorization_>>1));
4145            if ((progressFlag_&4)==0&&lastObj<thisObj+1.0e4&&
4146                largestPrimalError_<1.0e2) {
4147              // Just save costs
4148              // save extra copy of cost_
4149              int nTotal=numberRows_+numberColumns_;
4150              double * temp = new double [2*nTotal];
4151              memcpy(temp,cost_,nTotal*sizeof(double));
4152              memcpy(temp+nTotal,cost_,nTotal*sizeof(double));
4153              delete [] cost_;
4154              cost_=temp;
4155              objectiveWork_ = cost_;
4156              rowObjectiveWork_ = cost_+numberColumns_;
4157              progressFlag_ |= 4;
4158            } else {
4159              forceFactorization_=1;
4160#ifdef COIN_DEVELOP
4161              printf("Reducing factorization frequency - bad backwards\n");
4162#endif
4163#if 1
4164              unflagVariables = false;
4165              changeMade_++; // say something changed
4166              int nTotal = numberRows_+numberColumns_;
4167              CoinMemcpyN(saveStatus_,nTotal,status_);
4168              CoinMemcpyN(savedSolution_+numberColumns_ ,
4169                          numberRows_,rowActivityWork_);
4170              CoinMemcpyN(savedSolution_ ,
4171                          numberColumns_,columnActivityWork_);
4172              if ((progressFlag_&4)==0) {
4173                // save extra copy of cost_
4174                double * temp = new double [2*nTotal];
4175                memcpy(temp,cost_,nTotal*sizeof(double));
4176                memcpy(temp+nTotal,cost_,nTotal*sizeof(double));
4177                delete [] cost_;
4178                cost_=temp;
4179                objectiveWork_ = cost_;
4180                rowObjectiveWork_ = cost_+numberColumns_;
4181                progressFlag_ |= 4;
4182              } else {
4183                memcpy(cost_,cost_+nTotal,nTotal*sizeof(double));
4184              }
4185              // restore extra stuff
4186              int dummy;
4187              matrix_->generalExpanded(this,6,dummy);
4188              double pivotTolerance = factorization_->pivotTolerance();
4189              if(pivotTolerance<0.2)
4190                factorization_->pivotTolerance(0.2);
4191              else if(progress_.timesFlagged()>2)
4192                factorization_->pivotTolerance(CoinMin(pivotTolerance*1.1,0.99));
4193              if (alphaAccuracy_!=-1.0)
4194                alphaAccuracy_=-2.0;
4195              if (internalFactorize(1)) {
4196                CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
4197                CoinMemcpyN(savedSolution_+numberColumns_ ,
4198                            numberRows_,rowActivityWork_);
4199                CoinMemcpyN(savedSolution_ ,
4200                            numberColumns_,columnActivityWork_);
4201                // restore extra stuff
4202                int dummy;
4203                matrix_->generalExpanded(this,6,dummy);
4204                // debug
4205                int returnCode = internalFactorize(1);
4206                while (returnCode) {
4207                  // ouch
4208                  // switch off dense
4209                  int saveDense = factorization_->denseThreshold();
4210                  factorization_->setDenseThreshold(0);
4211                  // Go to safe
4212                  factorization_->pivotTolerance(0.99);
4213                  // make sure will do safe factorization
4214                  pivotVariable_[0]=-1;
4215                  returnCode=internalFactorize(2);
4216                  factorization_->setDenseThreshold(saveDense);
4217                }
4218              }
4219              resetFakeBounds(0);
4220              type = 2; // so will restore weights
4221              // get primal and dual solutions
4222              gutsOfSolution(givenDuals,NULL);
4223              if (numberPivots<2) {
4224                // need to reject something
4225                char x = isColumn(sequenceOut_) ? 'C' :'R';
4226                handler_->message(CLP_SIMPLEX_FLAG,messages_)
4227                  <<x<<sequenceWithin(sequenceOut_)
4228                  <<CoinMessageEol;
4229#ifdef COIN_DEVELOP
4230                printf("flag d\n");
4231#endif
4232                setFlagged(sequenceOut_);
4233                progress_.clearBadTimes();
4234                progress_.incrementTimesFlagged();
4235              }
4236              if (numberPivots<10)
4237                reallyBadProblems=true;
4238#ifdef COIN_DEVELOP
4239              printf("obj now %g\n",objectiveValue_);
4240#endif
4241              progress_.modifyObjective(objectiveValue_
4242                                        -bestPossibleImprovement_);
4243#endif
4244            }
4245          }
4246        } else {
4247          // in fast dual give up
4248#ifdef COIN_DEVELOP
4249          printf("In fast dual?\n");
4250#endif
4251          problemStatus_=3;
4252        }
4253      } else if (lastObj<thisObj-1.0e-5*CoinMax(fabs(thisObj),fabs(lastObj))-1.0e-3) {
4254        numberTimesOptimal_=0;
4255      }
4256    }
4257#endif
4258  }
4259  // Up tolerance if looks a bit odd
4260  if (numberIterations_>CoinMax(1000,numberRows_>>4)&&(specialOptions_&64)!=0) {
4261    if (sumPrimalInfeasibilities_&&sumPrimalInfeasibilities_<1.0e5) {
4262      int backIteration = progress_.lastIterationNumber(CLP_PROGRESS-1);
4263      if (backIteration>0&&numberIterations_-backIteration<9*CLP_PROGRESS) {
4264        if (factorization_->pivotTolerance()<0.9) {
4265          // up tolerance
4266          factorization_->pivotTolerance(CoinMin(factorization_->pivotTolerance()*1.05+0.02,0.91));
4267          //printf("tol now %g\n",factorization_->pivotTolerance());
4268          progress_.clearIterationNumbers();
4269        }
4270      }
4271    }
4272  }
4273  // Check if looping
4274  int loop;
4275  if (!givenDuals&&type!=2) 
4276    loop = progress_.looping();
4277  else
4278    loop=-1;
4279  if (progress_.reallyBadTimes()>10) {
4280    problemStatus_ = 10; // instead - try other algorithm
4281#if COIN_DEVELOP>2
4282    printf("returning at %d\n",__LINE__);
4283#endif
4284  }
4285  int situationChanged=0;
4286  if (loop>=0) {
4287    problemStatus_ = loop; //exit if in loop
4288    if (!problemStatus_) {
4289      // declaring victory
4290      numberPrimalInfeasibilities_ = 0;
4291      sumPrimalInfeasibilities_ = 0.0;
4292    } else {
4293      problemStatus_ = 10; // instead - try other algorithm
4294#if COIN_DEVELOP>2
4295      printf("returning at %d\n",__LINE__);
4296#endif
4297    }
4298    return;
4299  } else if (loop<-1) {
4300    // something may have changed
4301    gutsOfSolution(NULL,NULL);
4302    situationChanged=1;
4303  }
4304  // really for free variables in
4305  if((progressFlag_&2)!=0) {
4306    situationChanged=2;
4307  }
4308  progressFlag_ &= (~3); //reset progress flag
4309  if ((progressFlag_&4)!=0) {
4310    // save copy of cost_
4311    int nTotal = numberRows_+numberColumns_;
4312    memcpy(cost_+nTotal,cost_,nTotal*sizeof(double));
4313  }
4314  /*if (!numberIterations_&&sumDualInfeasibilities_)
4315    printf("OBJ %g sumPinf %g sumDinf %g\n",
4316           objectiveValue(),sumPrimalInfeasibilities_,
4317           sumDualInfeasibilities_);*/
4318  if (handler_->detail(CLP_SIMPLEX_STATUS,messages_)<100) {
4319    handler_->message(CLP_SIMPLEX_STATUS,messages_)
4320      <<numberIterations_<<objectiveValue();
4321    handler_->printing(sumPrimalInfeasibilities_>0.0)
4322      <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
4323    handler_->printing(sumDualInfeasibilities_>0.0)
4324      <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
4325    handler_->printing(numberDualInfeasibilitiesWithoutFree_
4326                       <numberDualInfeasibilities_)
4327                         <<numberDualInfeasibilitiesWithoutFree_;
4328    handler_->message()<<CoinMessageEol;
4329  }
4330#ifdef CLP_REPORT_PROGRESS
4331    if (ixxxxxx>=ixxyyyy-4&&ixxxxxx<=ixxyyyy) {
4332      char temp[20];
4333      sprintf(temp,"x_sol%d.out",ixxxxxx);
4334      FILE * fp = fopen(temp,"w");
4335      int nTotal=numberRows_+numberColumns_;
4336      for (int i=0;i<nTotal;i++)
4337        fprintf(fp,"%d %d %g %g %g %g %g\n",
4338                i,status_[i],lower_[i],solution_[i],upper_[i],cost_[i],dj_[i]);
4339      fclose(fp);
4340      if (ixxxxxx==ixxyyyy)
4341        exit(6);
4342    }
4343#endif
4344  realDualInfeasibilities=sumDualInfeasibilities_;
4345  double saveTolerance =dualTolerance_;
4346  // If we need to carry on cleaning variables
4347  if (!numberPrimalInfeasibilities_&&(specialOptions_&1024)!=0&&CLEAN_FIXED) {
4348    for (int iRow=0;iRow<numberRows_;iRow++) {
4349      int iPivot=pivotVariable_[iRow];
4350      if (!flagged(iPivot)&&pivoted(iPivot)) {
4351        // carry on
4352        numberPrimalInfeasibilities_=-1;
4353        sumOfRelaxedPrimalInfeasibilities_ = 1.0;
4354        sumPrimalInfeasibilities_ = 1.0;
4355        break;
4356      }
4357    }
4358  }
4359  /* If we are primal feasible and any dual infeasibilities are on
4360     free variables then it is better to go to primal */
4361  if (!numberPrimalInfeasibilities_&&!numberDualInfeasibilitiesWithoutFree_&&
4362      numberDualInfeasibilities_)
4363    problemStatus_=10;
4364  // dual bound coming in
4365  //double saveDualBound = dualBound_;
4366  bool needCleanFake=false;
4367  while (problemStatus_<=-3) {
4368    int cleanDuals=0;
4369    if (situationChanged!=0)
4370      cleanDuals=1;
4371    int numberChangedBounds=0;
4372    int doOriginalTolerance=0;
4373    if ( lastCleaned==numberIterations_)
4374      doOriginalTolerance=1;
4375    // check optimal
4376    // give code benefit of doubt
4377    if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
4378        sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
4379      // say optimal (with these bounds etc)
4380      numberDualInfeasibilities_ = 0;
4381      sumDualInfeasibilities_ = 0.0;
4382      numberPrimalInfeasibilities_ = 0;
4383      sumPrimalInfeasibilities_ = 0.0;
4384    }
4385    //if (dualFeasible()||problemStatus_==-4||(primalFeasible()&&!numberDualInfeasibilitiesWithoutFree_)) {
4386    if (dualFeasible()||problemStatus_==-4) {
4387      progress_.modifyObjective(objectiveValue_
4388                                -bestPossibleImprovement_);
4389#ifdef COIN_DEVELOP
4390      if (sumDualInfeasibilities_||bestPossibleImprovement_)
4391        printf("improve %g dualinf %g -> %g\n",
4392               bestPossibleImprovement_,sumDualInfeasibilities_,
4393               sumDualInfeasibilities_*dualBound_);
4394#endif
4395      // see if cutoff reached
4396      double limit = 0.0;
4397      getDblParam(ClpDualObjectiveLimit, limit);
4398#if 0
4399      if(fabs(limit)<1.0e30&&objectiveValue()*optimizationDirection_>
4400         limit+1.0e-7+1.0e-8*fabs(limit)&&!numberAtFakeBound()) {
4401        //looks infeasible on objective
4402        if (perturbation_==101) {
4403          cleanDuals=1;
4404          // Save costs
4405          int numberTotal = numberRows_+numberColumns_;
4406          double * saveCost = CoinCopyOfArray(cost_,numberTotal);
4407          // make sure fake bounds are back
4408          changeBounds(1,NULL,changeCost);
4409          createRim4(false);
4410          // make sure duals are current
4411          computeDuals(givenDuals);
4412          checkDualSolution();
4413          if(objectiveValue()*optimizationDirection_>
4414             limit+1.0e-7+1.0e-8*fabs(limit)&&!numberDualInfeasibilities_) {
4415            perturbation_=102; // stop any perturbations
4416            printf("cutoff test succeeded\n");
4417          } else {
4418            printf("cutoff test failed\n");
4419            // put back
4420            memcpy(cost_,saveCost,numberTotal*sizeof(double));
4421            // make sure duals are current
4422            computeDuals(givenDuals);
4423            checkDualSolution();
4424            progress_.modifyObjective(-COIN_DBL_MAX);
4425            problemStatus_=-1;
4426          }
4427          delete [] saveCost;
4428        }
4429      }
4430#endif
4431      if (primalFeasible()&&!givenDuals) {
4432        // may be optimal - or may be bounds are wrong
4433        handler_->message(CLP_DUAL_BOUNDS,messages_)
4434          <<dualBound_
4435          <<CoinMessageEol;
4436        // save solution in case unbounded
4437        double * saveColumnSolution = NULL;
4438        double * saveRowSolution = NULL;
4439        bool inCbc = (specialOptions_&(0x01000000|16384))!=0;
4440        if (!inCbc) {
4441          saveColumnSolution = CoinCopyOfArray(columnActivityWork_,numberColumns_);
4442          saveRowSolution = CoinCopyOfArray(rowActivityWork_,numberRows_);
4443        }
4444        numberChangedBounds=changeBounds(0,rowArray_[3],changeCost);
4445        if (numberChangedBounds<=0&&!numberDualInfeasibilities_) {
4446          //looks optimal - do we need to reset tolerance
4447          if (perturbation_==101) {
4448            perturbation_=102; // stop any perturbations
4449            cleanDuals=1;
4450            // make sure fake bounds are back
4451            //computeObjectiveValue();
4452            changeBounds(1,NULL,changeCost);
4453            //computeObjectiveValue();
4454            createRim4(false);
4455            // make sure duals are current
4456            computeDuals(givenDuals);
4457            checkDualSolution();
4458            progress_.modifyObjective(-COIN_DBL_MAX);
4459#define DUAL_TRY_FASTER
4460#ifdef DUAL_TRY_FASTER
4461            if (numberDualInfeasibilities_) {
4462#endif
4463              numberChanged_=1; // force something to happen
4464              lastCleaned=numberIterations_-1;
4465#ifdef DUAL_TRY_FASTER
4466            } else {
4467              //double value = objectiveValue_;
4468              computeObjectiveValue(true);
4469              //printf("old %g new %g\n",value,objectiveValue_);
4470              //numberChanged_=1;
4471            }
4472#endif
4473          }
4474          if (lastCleaned<numberIterations_&&numberTimesOptimal_<4&&
4475              (numberChanged_||(specialOptions_&4096)==0)) {
4476            doOriginalTolerance=2;
4477            numberTimesOptimal_++;
4478            changeMade_++; // say something changed
4479            if (numberTimesOptimal_==1) {
4480              dualTolerance_ = dblParam_[ClpDualTolerance];
4481            } else {
4482              if (numberTimesOptimal_==2) {
4483                // better to have small tolerance even if slower
4484                factorization_->zeroTolerance(CoinMin(factorization_->zeroTolerance(),1.0e-15));
4485              }
4486              dualTolerance_ = dblParam_[ClpDualTolerance];
4487              dualTolerance_ *= pow(2.0,numberTimesOptimal_-1);
4488            }
4489            cleanDuals=2; // If nothing changed optimal else primal
4490          } else {
4491            problemStatus_=0; // optimal
4492            if (lastCleaned<numberIterations_&&numberChanged_) {
4493              handler_->message(CLP_SIMPLEX_GIVINGUP,messages_)
4494                <<CoinMessageEol;
4495            }
4496          }
4497        } else {
4498          cleanDuals=1;
4499          if (doOriginalTolerance==1) {
4500            // check unbounded
4501            // find a variable with bad dj
4502            int iSequence;
4503            int iChosen=-1;
4504            if (!inCbc) {
4505              double largest = 100.0*primalTolerance_;
4506              for (iSequence=0;iSequence<numberRows_+numberColumns_;
4507                   iSequence++) {
4508                double djValue = dj_[iSequence];
4509                double originalLo = originalLower(iSequence);
4510                double originalUp = originalUpper(iSequence);
4511                if (fabs(djValue)>fabs(largest)) {
4512                  if (getStatus(iSequence)!=basic) {
4513                    if (djValue>0&&originalLo<-1.0e20) {
4514                      if (djValue>fabs(largest)) {
4515                        largest=djValue;
4516                        iChosen=iSequence;
4517                      }
4518                    } else if (djValue<0&&originalUp>1.0e20) {
4519                      if (-djValue>fabs(largest)) {
4520                        largest=djValue;
4521                        iChosen=iSequence;
4522                      }
4523                    } 
4524                  }
4525                }
4526              }
4527            }
4528            if (iChosen>=0) {
4529              int iSave=sequenceIn_;
4530              sequenceIn_=iChosen;
4531              unpack(rowArray_[1]);
4532              sequenceIn_ = iSave;
4533              // if dual infeasibilities then must be free vector so add in dual
4534              if (numberDualInfeasibilities_) {
4535                if (fabs(changeCost)>1.0e-5)
4536                  printf("Odd free/unbounded combo\n");
4537                changeCost += cost_[iChosen];
4538              }
4539              problemStatus_ = checkUnbounded(rowArray_[1],rowArray_[0],
4540                                              changeCost);
4541              rowArray_[1]->clear();
4542            } else {
4543              problemStatus_=-3;
4544            }
4545            if (problemStatus_==2&&perturbation_==101) {
4546              perturbation_=102; // stop any perturbations
4547              cleanDuals=1;
4548              createRim4(false);
4549              progress_.modifyObjective(-COIN_DBL_MAX);
4550              problemStatus_=-1;
4551            }
4552            if (problemStatus_==2) {
4553              // it is unbounded - restore solution
4554              // but first add in changes to non-basic
4555              int iColumn;
4556              double * original = columnArray_[0]->denseVector();
4557              for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4558                if(getColumnStatus(iColumn)!= basic)
4559                  ray_[iColumn] += 
4560                    saveColumnSolution[iColumn]-original[iColumn];
4561                columnActivityWork_[iColumn] = original[iColumn];
4562              }
4563              CoinMemcpyN(saveRowSolution,numberRows_,
4564                                rowActivityWork_);
4565            }
4566          } else {
4567            doOriginalTolerance=2;
4568            rowArray_[0]->clear();
4569          }
4570        }
4571        delete [] saveColumnSolution;
4572        delete [] saveRowSolution;
4573      }
4574      if (problemStatus_==-4||problemStatus_==-5) {
4575        // may be infeasible - or may be bounds are wrong
4576        numberChangedBounds=changeBounds(0,NULL,changeCost);
4577        needCleanFake=true;
4578        /* Should this be here as makes no difference to being feasible.
4579           But seems to make a difference to run times. */
4580        if (perturbation_==101&&0) {
4581          perturbation_=102; // stop any perturbations
4582          cleanDuals=1;
4583          numberChangedBounds=1;
4584          // make sure fake bounds are back
4585          changeBounds(1,NULL,changeCost);
4586          needCleanFake=true;
4587          createRim4(false);
4588          progress_.modifyObjective(-COIN_DBL_MAX);
4589        }
4590        if ((numberChangedBounds<=0||dualBound_>1.0e20||
4591            (largestPrimalError_>1.0&&dualBound_>1.0e17))&&
4592            (numberPivots<4||sumPrimalInfeasibilities_>1.0e-6)) {
4593          problemStatus_=1; // infeasible
4594          if (perturbation_==101) {
4595            perturbation_=102; // stop any perturbations
4596            //cleanDuals=1;
4597            //numberChangedBounds=1;
4598            //createRim4(false);
4599          }
4600        } else {
4601          problemStatus_=-1; //iterate
4602          cleanDuals=1;
4603          if (numberChangedBounds<=0)
4604            doOriginalTolerance=2;
4605          // and delete ray which has been created
4606          delete [] ray_;
4607          ray_ = NULL;
4608        }
4609
4610      }
4611    } else {
4612      cleanDuals=1;
4613    }
4614    if (problemStatus_<0) {
4615      if (doOriginalTolerance==2) {
4616        // put back original tolerance
4617        lastCleaned=numberIterations_;
4618        numberChanged_ =0; // Number of variables with changed costs
4619        handler_->message(CLP_DUAL_ORIGINAL,messages_)
4620          <<CoinMessageEol;
4621        perturbation_=102; // stop any perturbations
4622#if 0
4623        double * xcost = new double[numberRows_+numberColumns_];
4624        double * xlower = new double[numberRows_+numberColumns_];
4625        double * xupper = new double[numberRows_+numberColumns_];
4626        double * xdj = new double[numberRows_+numberColumns_];
4627        double * xsolution = new double[numberRows_+numberColumns_];
4628 CoinMemcpyN(cost_,(numberRows_+numberColumns_),xcost);
4629 CoinMemcpyN(lower_,(numberRows_+numberColumns_),xlower);
4630 CoinMemcpyN(upper_,(numberRows_+numberColumns_),xupper);
4631 CoinMemcpyN(dj_,(numberRows_+numberColumns_),xdj);
4632 CoinMemcpyN(solution_,(numberRows_+numberColumns_),xsolution);
4633#endif
4634        createRim4(false);
4635        progress_.modifyObjective(-COIN_DBL_MAX);
4636        // make sure duals are current
4637        computeDuals(givenDuals);
4638        checkDualSolution();
4639#if 0
4640        int i;
4641        for (i=0;i<numberRows_+numberColumns_;i++) {
4642          if (cost_[i]!=xcost[i])
4643            printf("** %d old cost %g new %g sol %g\n",
4644                   i,xcost[i],cost_[i],solution_[i]);
4645          if (lower_[i]!=xlower[i])
4646            printf("** %d old lower %g new %g sol %g\n",
4647                   i,xlower[i],lower_[i],solution_[i]);
4648          if (upper_[i]!=xupper[i])
4649            printf("** %d old upper %g new %g sol %g\n",
4650                   i,xupper[i],upper_[i],solution_[i]);
4651          if (dj_[i]!=xdj[i])
4652            printf("** %d old dj %g new %g sol %g\n",
4653                   i,xdj[i],dj_[i],solution_[i]);
4654          if (solution_[i]!=xsolution[i])
4655            printf("** %d old solution %g new %g sol %g\n",
4656                   i,xsolution[i],solution_[i],solution_[i]);
4657        }
4658        //delete [] xcost;
4659        //delete [] xupper;
4660        //delete [] xlower;
4661        //delete [] xdj;
4662        //delete [] xsolution;
4663#endif
4664        // put back bounds as they were if was optimal
4665        if (doOriginalTolerance==2&&cleanDuals!=2) {
4666          changeMade_++; // say something changed
4667          /* We may have already changed some bounds in this function
4668             so save numberFake_ and add in.
4669
4670             Worst that can happen is that we waste a bit of time  - but it must be finite.
4671          */
4672          //int saveNumberFake = numberFake_;
4673          //resetFakeBounds(-1);
4674          changeBounds(3,NULL,changeCost);
4675          needCleanFake=true;
4676          //numberFake_ += saveNumberFake;
4677          //resetFakeBounds(-1);
4678          cleanDuals=2;
4679          //cleanDuals=1;
4680        }
4681#if 0
4682        //int i;
4683        for (i=0;i<numberRows_+numberColumns_;i++) {
4684          if (cost_[i]!=xcost[i])
4685            printf("** %d old cost %g new %g sol %g\n",
4686                   i,xcost[i],cost_[i],solution_[i]);
4687          if (lower_[i]!=xlower[i])
4688            printf("** %d old lower %g new %g sol %g\n",
4689                   i,xlower[i],lower_[i],solution_[i]);
4690          if (upper_[i]!=xupper[i])
4691            printf("** %d old upper %g new %g sol %g\n",
4692                   i,xupper[i],upper_[i],solution_[i]);
4693          if (dj_[i]!=xdj[i])
4694            printf("** %d old dj %g new %g sol %g\n",
4695                   i,xdj[i],dj_[i],solution_[i]);
4696          if (solution_[i]!=xsolution[i])
4697            printf("** %d old solution %g new %g sol %g\n",
4698                   i,xsolution[i],solution_[i],solution_[i]);
4699        }
4700        delete [] xcost;
4701        delete [] xupper;
4702        delete [] xlower;
4703        delete [] xdj;
4704        delete [] xsolution;
4705#endif
4706      }
4707      if (cleanDuals==1||(cleanDuals==2&&!numberDualInfeasibilities_)) {
4708        // make sure dual feasible
4709        // look at all rows and columns
4710        rowArray_[0]->clear();
4711        columnArray_[0]->clear();
4712        double objectiveChange=0.0;
4713#if 0
4714        double * xcost = new double[numberRows_+numberColumns_];
4715        double * xlower = new double[numberRows_+numberColumns_];
4716        double * xupper = new double[numberRows_+numberColumns_];
4717        double * xdj = new double[numberRows_+numberColumns_];
4718        double * xsolution = new double[numberRows_+numberColumns_];
4719 CoinMemcpyN(cost_,(numberRows_+numberColumns_),xcost);
4720 CoinMemcpyN(lower_,(numberRows_+numberColumns_),xlower);
4721 CoinMemcpyN(upper_,(numberRows_+numberColumns_),xupper);
4722 CoinMemcpyN(dj_,(numberRows_+numberColumns_),xdj);
4723 CoinMemcpyN(solution_,(numberRows_+numberColumns_),xsolution);
4724#endif
4725        if (givenDuals)
4726          dualTolerance_=1.0e50;
4727        updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[1],
4728          0.0,objectiveChange,true);
4729        dualTolerance_=saveTolerance;
4730#if 0
4731        int i;
4732        for (i=0;i<numberRows_+numberColumns_;i++) {
4733          if (cost_[i]!=xcost[i])
4734            printf("** %d old cost %g new %g sol %g\n",
4735                   i,xcost[i],cost_[i],solution_[i]);
4736          if (lower_[i]!=xlower[i])
4737            printf("** %d old lower %g new %g sol %g\n",
4738                   i,xlower[i],lower_[i],solution_[i]);
4739          if (upper_[i]!=xupper[i])
4740            printf("** %d old upper %g new %g sol %g\n",
4741                   i,xupper[i],upper_[i],solution_[i]);
4742          if (dj_[i]!=xdj[i])
4743            printf("** %d old dj %g new %g sol %g\n",
4744                   i,xdj[i],dj_[i],solution_[i]);
4745          if (solution_[i]!=xsolution[i])
4746            printf("** %d old solution %g new %g sol %g\n",
4747                   i,xsolution[i],solution_[i],solution_[i]);
4748        }
4749        delete [] xcost;
4750        delete [] xupper;
4751        delete [] xlower;
4752        delete [] xdj;
4753        delete [] xsolution;
4754#endif
4755        // for now - recompute all
4756        gutsOfSolution(NULL,NULL);
4757        if (givenDuals)
4758          dualTolerance_=1.0e50;
4759        updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[1],
4760          0.0,objectiveChange,true);
4761        dualTolerance_=saveTolerance;
4762        //assert(numberDualInfeasibilitiesWithoutFree_==0);
4763        if (numberDualInfeasibilities_) {
4764          if (numberPrimalInfeasibilities_||numberPivots) {
4765            problemStatus_=-1; // carry on as normal
4766          } else {
4767            problemStatus_=10; // try primal
4768#if COIN_DEVELOP>1
4769            printf("returning at %d\n",__LINE__);
4770#endif
4771          }
4772        } else if (situationChanged==2) {
4773          problemStatus_=-1; // carry on as normal
4774        }
4775        situationChanged=0;
4776      } else {
4777        // iterate
4778        if (cleanDuals!=2) {
4779          problemStatus_=-1;
4780        } else {
4781          problemStatus_ = 10; // try primal
4782#if COIN_DEVELOP>2
4783          printf("returning at %d\n",__LINE__);
4784#endif
4785        }
4786      }
4787    }
4788  }
4789  // unflag all variables (we may want to wait a bit?)
4790  if ((tentativeStatus!=-2&&tentativeStatus!=-1)&&unflagVariables) {
4791    int iRow;
4792    int numberFlagged=0;
4793    for (iRow=0;iRow<numberRows_;iRow++) {
4794      int iPivot=pivotVariable_[iRow];
4795      if (flagged(iPivot)) {
4796        numberFlagged++;
4797        clearFlagged(iPivot);
4798      }
4799    }
4800#ifdef COIN_DEVELOP
4801    if (numberFlagged) {
4802      printf("unflagging %d variables - tentativeStatus %d probStat %d ninf %d nopt %d\n",numberFlagged,tentativeStatus,
4803             problemStatus_,numberPrimalInfeasibilities_,
4804             numberTimesOptimal_);
4805    }
4806#endif
4807    unflagVariables = numberFlagged>0;
4808    if (numberFlagged&&!numberPivots) {
4809      /* looks like trouble as we have not done any iterations.
4810         Try changing pivot tolerance then give it a few goes and give up */
4811      if (factorization_->pivotTolerance()<0.9) {
4812        factorization_->pivotTolerance(0.99);
4813        problemStatus_=-1;
4814      } else if (numberTimesOptimal_<3) {
4815        numberTimesOptimal_++;
4816        problemStatus_=-1;
4817      } else {
4818        unflagVariables=false;
4819        //secondaryStatus_ = 1; // and say probably infeasible
4820        // try primal
4821        problemStatus_=10;
4822#if COIN_DEVELOP>1
4823        printf("returning at %d\n",__LINE__);
4824#endif
4825      }
4826    }
4827  }
4828  if (problemStatus_<0) {
4829    if (needCleanFake) {
4830      double dummyChangeCost=0.0;
4831      changeBounds(3,NULL,dummyChangeCost);
4832    }
4833#if 0
4834    if (objectiveValue_<lastObjectiveValue_-1.0e-8*
4835        CoinMax(fabs(objectivevalue_),fabs(lastObjectiveValue_))) {
4836    } else {
4837      lastObjectiveValue_=objectiveValue_;
4838    }
4839#endif
4840    if (type==0||type==1) {
4841      if (!type) {
4842        // create save arrays
4843        delete [] saveStatus_;
4844        delete [] savedSolution_;
4845        saveStatus_ = new unsigned char [numberRows_+numberColumns_];
4846        savedSolution_ = new double [numberRows_+numberColumns_];
4847      }
4848      // save arrays
4849      CoinMemcpyN(status_,numberColumns_+numberRows_,saveStatus_);
4850      CoinMemcpyN(rowActivityWork_,
4851                  numberRows_,savedSolution_+numberColumns_);
4852      CoinMemcpyN(columnActivityWork_,numberColumns_,savedSolution_);
4853      // save extra stuff
4854      int dummy;
4855      matrix_->generalExpanded(this,5,dummy);
4856    }
4857    if (weightsSaved) {
4858      // restore weights (if saved) - also recompute infeasibility list
4859      if (!reallyBadProblems&&(largestPrimalError_<100.0||numberPivots>10)) {
4860        if (tentativeStatus>-3) 
4861          dualRowPivot_->saveWeights(this,(type <2) ? 2 : 4);
4862        else
4863          dualRowPivot_->saveWeights(this,3);
4864      } else {
4865        // reset weights or scale back
4866        dualRowPivot_->saveWeights(this,6);
4867      }
4868    }
4869  }
4870  // see if cutoff reached
4871  double limit = 0.0;
4872  getDblParam(ClpDualObjectiveLimit, limit);
4873#if 0
4874  if(fabs(limit)<1.0e30&&objectiveValue()*optimizationDirection_>
4875     limit+100.0) {
4876    printf("lim %g obj %g %g - wo perturb %g sum dual %g\n",
4877           limit,objectiveValue_,objectiveValue(),computeInternalObjectiveValue(),sumDualInfeasibilities_);
4878  }
4879#endif
4880  if(fabs(limit)<1.0e30&&objectiveValue()*optimizationDirection_>
4881     limit&&!numberAtFakeBound()) {
4882    bool looksInfeasible = !numberDualInfeasibilities_;
4883    if (objectiveValue()*optimizationDirection_>limit+fabs(0.1*limit)+1.0e2*sumDualInfeasibilities_+1.0e4&&
4884        sumDualInfeasibilities_<largestDualError_&&numberIterations_>0.5*numberRows_+1000)
4885      looksInfeasible=true;
4886    if (looksInfeasible) {
4887      // Even if not perturbed internal costs may have changed
4888      // be careful
4889      if (true||numberIterations_) {
4890        if(computeInternalObjectiveValue()>limit) {
4891          problemStatus_=1;
4892          secondaryStatus_ = 1; // and say was on cutoff
4893        }
4894      } else {
4895        problemStatus_=1;
4896        secondaryStatus_ = 1; // and say was on cutoff
4897      }
4898    }
4899  }
4900  // If we are in trouble and in branch and bound give up
4901  if ((specialOptions_&1024)!=0) {
4902    int looksBad=0;
4903    if (largestPrimalError_*largestDualError_>1.0e2) {
4904      looksBad=1;
4905    } else if (largestPrimalError_>1.0e-2
4906        &&objectiveValue_>CoinMin(1.0e15,1.0e3*limit)) {
4907      looksBad=2;
4908    }
4909    if (looksBad) {
4910      if (factorization_->pivotTolerance()<0.9) {
4911        // up tolerance
4912        factorization_->pivotTolerance(CoinMin(factorization_->pivotTolerance()*1.05+0.02,0.91));
4913      } else if (numberIterations_>10000) {
4914        if (handler_->logLevel()>2)
4915          printf("bad dual - saying infeasible %d\n",looksBad);
4916        problemStatus_=1;
4917        secondaryStatus_ = 1; // and say was on cutoff
4918      } else if (largestPrimalError_>1.0e5) {
4919        {
4920          int iBigB=-1;
4921          double bigB=0.0;
4922          int iBigN=-1;
4923          double bigN=0.0;
4924          for (int i=0;i<numberRows_+numberColumns_;i++) {
4925            double value = fabs(solution_[i]);
4926            if (getStatus(i)==basic) {
4927              if (value>bigB) {
4928                bigB= value;
4929                iBigB=i;
4930              }
4931            } else {
4932              if (value>bigN) {
4933                bigN= value;
4934                iBigN=i;
4935              }
4936            }
4937          }
4938#ifdef CLP_INVESTIGATE
4939          if (bigB>1.0e8||bigN>1.0e8) {
4940            if (handler_->logLevel()>0)
4941              printf("it %d - basic %d %g, nonbasic %d %g\n",
4942                     numberIterations_,iBigB,bigB,iBigN,bigN);
4943          }
4944#endif
4945        }
4946#if COIN_DEVELOP!=2
4947        if (handler_->logLevel()>2)
4948#endif
4949          printf("bad dual - going to primal %d %g\n",looksBad,largestPrimalError_);
4950        allSlackBasis(true);
4951        problemStatus_=10;
4952      }
4953    }
4954  }
4955  if (problemStatus_<0&&!changeMade_) {
4956    problemStatus_=4; // unknown
4957  }
4958  lastGoodIteration_ = numberIterations_;
4959  if (numberIterations_>lastBadIteration_+100)
4960    moreSpecialOptions_ &= ~16; // clear check accuracy flag
4961  if (problemStatus_<0) {
4962    sumDualInfeasibilities_=realDualInfeasibilities; // back to say be careful
4963    if (sumDualInfeasibilities_)
4964      numberDualInfeasibilities_=1;
4965  }
4966#ifdef CLP_REPORT_PROGRESS
4967  if (ixxxxxx>ixxyyyy-3) {
4968    printf("objectiveValue_ %g\n",objectiveValue_);
4969    handler_->setLogLevel(63);
4970    int nTotal = numberColumns_+numberRows_;
4971    double newObj=0.0;
4972    for (int i=0;i<nTotal;i++) {
4973      if (solution_[i])
4974        newObj += solution_[i]*cost_[i];
4975    }
4976    printf("xxx obj %g\n",newObj);
4977    // for now - recompute all
4978    gutsOfSolution(NULL,NULL);
4979    newObj=0.0;
4980    for (int i=0;i<nTotal;i++) {
4981      if (solution_[i])
4982        newObj += solution_[i]*cost_[i];
4983    }
4984    printf("yyy obj %g %g\n",newObj,objectiveValue_);
4985    progress_.modifyObjective(objectiveValue_
4986                              -bestPossibleImprovement_);
4987  }
4988#endif
4989#if 1
4990  double thisObj = progress_.lastObjective(0);
4991  double lastObj = progress_.lastObjective(1);
4992  if (lastObj>thisObj+1.0e-4*CoinMax(fabs(thisObj),fabs(lastObj))+1.0e-4
4993      &&givenDuals==NULL&&firstFree_<0) {
4994    int maxFactor = factorization_->maximumPivots();
4995    if (maxFactor>10) {
4996      if (forceFactorization_<0)
4997        forceFactorization_= maxFactor;
4998      forceFactorization_ = CoinMax(1,(forceFactorization_>>1));
4999      //printf("Reducing factorization frequency\n");
5000    } 
5001  }
5002#endif
5003  // Allow matrices to be sorted etc
5004  int fake=-999; // signal sort
5005  matrix_->correctSequence(this,fake,fake);
5006  if (alphaAccuracy_>0.0)
5007      alphaAccuracy_=1.0;
5008}
5009/* While updateDualsInDual sees what effect is of flip
5010   this does actual flipping.
5011   If change >0.0 then value in array >0.0 => from lower to upper
5012*/
5013void 
5014ClpSimplexDual::flipBounds(CoinIndexedVector * rowArray,
5015                  CoinIndexedVector * columnArray,
5016                  double change)
5017{
5018  int number;
5019  int * which;
5020 
5021  int iSection;
5022
5023  for (iSection=0;iSection<2;iSection++) {
5024    int i;
5025    double * solution = solutionRegion(iSection);
5026    double * lower = lowerRegion(iSection);
5027    double * upper = upperRegion(iSection);
5028    int addSequence;
5029    if (!iSection) {
5030      number = rowArray->getNumElements();
5031      which = rowArray->getIndices();
5032      addSequence = numberColumns_;
5033    } else {
5034      number = columnArray->getNumElements();
5035      which = columnArray->getIndices();
5036      addSequence = 0;
5037    }
5038   
5039    for (i=0;i<number;i++) {
5040      int iSequence = which[i];
5041      Status status = getStatus(iSequence+addSequence);
5042
5043      switch(status) {
5044
5045      case basic:
5046      case isFree:
5047      case superBasic:
5048      case ClpSimplex::isFixed:
5049        break;
5050      case atUpperBound:
5051        // to lower bound
5052        setStatus(iSequence+addSequence,atLowerBound);
5053        solution[iSequence] = lower[iSequence];
5054        break;
5055      case atLowerBound:
5056        // to upper bound
5057        setStatus(iSequence+addSequence,atUpperBound);
5058        solution[iSequence] = upper[iSequence];
5059        break;
5060      }
5061    }
5062  }
5063  rowArray->setNumElements(0);
5064  columnArray->setNumElements(0);
5065}
5066// Restores bound to original bound
5067void 
5068ClpSimplexDual::originalBound( int iSequence)
5069{
5070  if (getFakeBound(iSequence)!=noFake) {
5071    numberFake_--;;
5072    setFakeBound(iSequence,noFake);
5073#ifdef CLP_AUXILIARY_MODEL
5074    if (auxiliaryModel_) {
5075      // just copy back
5076      lower_[iSequence]=auxiliaryModel_->lowerRegion()[iSequence+numberRows_+numberColumns_];
5077      upper_[iSequence]=auxiliaryModel_->upperRegion()[iSequence+numberRows_+numberColumns_];
5078      return;
5079    }
5080#endif
5081    if (iSequence>=numberColumns_) {
5082      // rows
5083      int iRow = iSequence-numberColumns_;
5084      rowLowerWork_[iRow]=rowLower_[iRow];
5085      rowUpperWork_[iRow]=rowUpper_[iRow];
5086      if (rowScale_) {
5087        if (rowLowerWork_[iRow]>-1.0e50)
5088          rowLowerWork_[iRow] *= rowScale_[iRow]*rhsScale_;
5089        if (rowUpperWork_[iRow]<1.0e50)
5090          rowUpperWork_[iRow] *= rowScale_[iRow]*rhsScale_;
5091      } else if (rhsScale_!=1.0) {
5092        if (rowLowerWork_[iRow]>-1.0e50)
5093          rowLowerWork_[iRow] *= rhsScale_;
5094        if (rowUpperWork_[iRow]<1.0e50)
5095          rowUpperWork_[iRow] *= rhsScale_;
5096      }
5097    } else {
5098      // columns
5099      columnLowerWork_[iSequence]=columnLower_[iSequence];
5100      columnUpperWork_[iSequence]=columnUpper_[iSequence];
5101      if (rowScale_) {
5102        double multiplier = 1.0/columnScale_[iSequence];
5103        if (columnLowerWork_[iSequence]>-1.0e50)
5104          columnLowerWork_[iSequence] *= multiplier*rhsScale_;
5105        if (columnUpperWork_[iSequence]<1.0e50)
5106          columnUpperWork_[iSequence] *= multiplier*rhsScale_;
5107      } else if (rhsScale_!=1.0) {
5108        if (columnLowerWork_[iSequence]>-1.0e50)
5109          columnLowerWork_[iSequence] *= rhsScale_;
5110        if (columnUpperWork_[iSequence]<1.0e50)
5111          columnUpperWork_[iSequence] *= rhsScale_;
5112      }
5113    }
5114  }
5115}
5116/* As changeBounds but just changes new bounds for a single variable.
5117   Returns true if change */
5118bool 
5119ClpSimplexDual::changeBound( int iSequence)
5120{
5121  // old values
5122  double oldLower=lower_[iSequence];
5123  double oldUpper=upper_[iSequence];
5124  double value=solution_[iSequence];
5125  bool modified=false;
5126  originalBound(iSequence);
5127  // original values
5128  double lowerValue=lower_[iSequence];
5129  double upperValue=upper_[iSequence];
5130  // back to altered values
5131  lower_[iSequence] = oldLower;
5132  upper_[iSequence] = oldUpper;
5133  assert (getFakeBound(iSequence)==noFake);
5134  //if (getFakeBound(iSequence)!=noFake)
5135  //numberFake_--;;
5136  if (value==oldLower) {
5137    if (upperValue > oldLower + dualBound_) {
5138      upper_[iSequence]=oldLower+dualBound_;
5139      setFakeBound(iSequence,upperFake);
5140      modified=true;
5141      numberFake_++;
5142    }
5143  } else if (value==oldUpper) {
5144    if (lowerValue < oldUpper - dualBound_) {
5145      lower_[iSequence]=oldUpper-dualBound_;
5146      setFakeBound(iSequence,lowerFake);
5147      modified=true;
5148      numberFake_++;
5149    }
5150  } else {
5151    assert(value==oldLower||value==oldUpper);
5152  }
5153  return modified;
5154}
5155// Perturbs problem
5156int 
5157ClpSimplexDual::perturb()
5158{
5159  if (perturbation_>100)
5160    return 0; //perturbed already
5161  if (perturbation_==100)
5162    perturbation_=50; // treat as normal
5163  int savePerturbation = perturbation_;
5164  bool modifyRowCosts=false;
5165  // dual perturbation
5166  double perturbation=1.0e-20;
5167  // maximum fraction of cost to perturb
5168  double maximumFraction = 1.0e-5;
5169  double constantPerturbation = 100.0*dualTolerance_;
5170  int maxLength=0;
5171  int minLength=numberRows_;
5172  double averageCost = 0.0;
5173#if 0
5174  // look at element range
5175  double smallestNegative;
5176  double largestNegative;
5177  double smallestPositive;
5178  double largestPositive;
5179  matrix_->rangeOfElements(smallestNegative, largestNegative,
5180                           smallestPositive, largestPositive);
5181  smallestPositive = CoinMin(fabs(smallestNegative),smallestPositive);
5182  largestPositive = CoinMax(fabs(largestNegative),largestPositive);
5183  double elementRatio = largestPositive/smallestPositive;
5184#endif
5185  int numberNonZero=0;
5186  if (!numberIterations_&&perturbation_>=50) {
5187    // See if we need to perturb
5188    double * sort = new double[numberColumns_];
5189    // Use objective BEFORE scaling
5190    const double * obj = objective();
5191    int i;
5192    for (i=0;i<numberColumns_;i++) {
5193      double value = fabs(obj[i]);
5194      sort[i]=value;
5195      averageCost += value;
5196      if (value)
5197        numberNonZero++;
5198    }
5199    if (numberNonZero)
5200      averageCost /= static_cast<double> (numberNonZero);
5201    else
5202      averageCost = 1.0;
5203    std::sort(sort,sort+numberColumns_);
5204    int number=1;
5205    double last = sort[0];
5206    for (i=1;i<numberColumns_;i++) {
5207      if (last!=sort[i])
5208        number++;
5209      last=sort[i];
5210    }
5211    delete [] sort;
5212    if (!numberNonZero)
5213      return 1; // safer to use primal
5214#if 0
5215    printf("nnz %d percent %d",number,(number*100)/numberColumns_);
5216    if (number*4>numberColumns_)
5217      printf(" - Would not perturb\n");
5218    else
5219      printf(" - Would perturb\n");
5220    //exit(0);
5221#endif
5222    //printf("ratio number diff costs %g, element ratio %g\n",((double)number)/((double) numberColumns_),
5223    //                                                                elementRatio);
5224    //number=0;
5225    //if (number*4>numberColumns_||elementRatio>1.0e12) {
5226    if (number*4>numberColumns_) {
5227      perturbation_=100;
5228      return 0; // good enough
5229    }
5230  }
5231  int iColumn;
5232  const int * columnLength = matrix_->getVectorLengths();
5233  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
5234    if (columnLowerWork_[iColumn]<columnUpperWork_[iColumn]) {
5235      int length = columnLength[iColumn];
5236      if (length>2) {
5237        maxLength = CoinMax(maxLength,length);
5238        minLength = CoinMin(minLength,length);
5239      }
5240    }
5241  }
5242  // If > 70 then do rows
5243  if (perturbation_>=70) {
5244    modifyRowCosts=true;
5245    perturbation_ -= 20;
5246    printf("Row costs modified, ");
5247  }
5248  bool uniformChange=false;
5249  bool inCbcOrOther = (specialOptions_&0x03000000)!=0;
5250  if (perturbation_>50) {
5251    // Experiment
5252    // maximumFraction could be 1.0e-10 to 1.0
5253    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};
5254    int whichOne = perturbation_-51;
5255    //if (inCbcOrOther&&whichOne>0)
5256    //whichOne--;
5257    maximumFraction = m[CoinMin(whichOne,10)];
5258  } else if (inCbcOrOther) {
5259    //maximumFraction = 1.0e-6;
5260  }
5261  int iRow;
5262  double smallestNonZero=1.0e100;
5263  numberNonZero=0;
5264  if (perturbation_>=50) {
5265    perturbation = 1.0e-8;
5266    bool allSame=true;
5267    double lastValue=0.0;
5268    for (iRow=0;iRow<numberRows_;iRow++) {
5269      double lo = rowLowerWork_[iRow];
5270      double up = rowUpperWork_[iRow];
5271      if (lo<up) {
5272        double value = fabs(rowObjectiveWork_[iRow]);
5273        perturbation = CoinMax(perturbation,value);
5274        if (value) {
5275          modifyRowCosts=true;
5276          smallestNonZero = CoinMin(smallestNonZero,value);
5277        }
5278      } 
5279      if (lo&&lo>-1.0e10) {
5280        numberNonZero++;
5281        lo=fabs(lo);
5282        if (!lastValue) 
5283          lastValue=lo;
5284        else if (fabs(lo-lastValue)>1.0e-7)
5285          allSame=false;
5286      }
5287      if (up&&up<1.0e10) {
5288        numberNonZero++;
5289        up=fabs(up);
5290        if (!lastValue) 
5291          lastValue=up;
5292        else if (fabs(up-lastValue)>1.0e-7)
5293          allSame=false;
5294      }
5295    }
5296    double lastValue2=0.0;
5297    for (iColumn=0;iColumn<numberColumns_;iColumn++) { 
5298      double lo = columnLowerWork_[iColumn];
5299      double up = columnUpperWork_[iColumn];
5300      if (lo<up) {
5301        double value = 
5302          fabs(objectiveWork_[iColumn]);
5303        perturbation = CoinMax(perturbation,value);
5304        if (value) {
5305          smallestNonZero = CoinMin(smallestNonZero,value);
5306        }
5307      }
5308      if (lo&&lo>-1.0e10) {
5309        //numberNonZero++;
5310        lo=fabs(lo);
5311        if (!lastValue2) 
5312          lastValue2=lo;
5313        else if (fabs(lo-lastValue2)>1.0e-7)
5314          allSame=false;
5315      }
5316      if (up&&up<1.0e10) {
5317        //numberNonZero++;
5318        up=fabs(up);
5319        if (!lastValue2) 
5320          lastValue2=up;
5321        else if (fabs(up-lastValue2)>1.0e-7)
5322          allSame=false;
5323      }
5324    }
5325    if (allSame) {
5326      // Check elements
5327      double smallestNegative;
5328      double largestNegative;
5329      double smallestPositive;
5330      double largestPositive;
5331      matrix_->rangeOfElements(smallestNegative,largestNegative,
5332                               smallestPositive,largestPositive);
5333      if (smallestNegative==largestNegative&&
5334          smallestPositive==largestPositive) {
5335        // Really hit perturbation
5336        double adjust = CoinMin(100.0*maximumFraction,1.0e-3*CoinMax(lastValue,lastValue2));
5337        maximumFraction = CoinMax(adjust,maximumFraction);
5338      }
5339    }
5340    perturbation = CoinMin(perturbation,smallestNonZero/maximumFraction);
5341  } else {
5342    // user is in charge
5343    maximumFraction = 1.0e-1;
5344    // but some experiments
5345    if (perturbation_<=-900) {
5346      modifyRowCosts=true;
5347      perturbation_ += 1000;
5348      printf("Row costs modified, ");
5349    }
5350    if (perturbation_<=-10) {
5351      perturbation_ += 10; 
5352      maximumFraction = 1.0;
5353      if ((-perturbation_)%100>=10) {
5354        uniformChange=true;
5355        perturbation_+=20;
5356      }
5357      while (perturbation_<-10) {
5358        perturbation_ += 100;
5359        maximumFraction *= 1.0e-1;
5360      }
5361    }
5362    perturbation = pow(10.0,perturbation_);
5363  }
5364  double largestZero=0.0;
5365  double largest=0.0;
5366  double largestPerCent=0.0;
5367  // modify costs
5368  bool printOut=(handler_->logLevel()==63);
5369  printOut=false;
5370  //assert (!modifyRowCosts);
5371  modifyRowCosts=false;
5372  if (modifyRowCosts) {
5373    for (iRow=0;iRow<numberRows_;iRow++) {
5374      if (rowLowerWork_[iRow]<rowUpperWork_[iRow]) {
5375        double value = perturbation;
5376        double currentValue = rowObjectiveWork_[iRow];
5377        value = CoinMin(value,maximumFraction*(fabs(currentValue)+1.0e-1*perturbation+1.0e-3));
5378        if (rowLowerWork_[iRow]>-largeValue_) {
5379          if (fabs(rowLowerWork_[iRow])<fabs(rowUpperWork_[iRow])) 
5380            value *= randomNumberGenerator_.randomDouble();
5381          else
5382            value *= -randomNumberGenerator_.randomDouble();
5383        } else if (rowUpperWork_[iRow]<largeValue_) {
5384          value *= -randomNumberGenerator_.randomDouble();
5385        } else {
5386          value=0.0;
5387        }
5388        if (currentValue) {
5389          largest = CoinMax(largest,fabs(value));
5390          if (fabs(value)>fabs(currentValue)*largestPerCent)
5391            largestPerCent=fabs(value/currentValue);
5392        } else {
5393          largestZero=CoinMax(largestZero,fabs(value));
5394        }
5395        if (printOut)
5396          printf("row %d cost %g change %g\n",iRow,rowObjectiveWork_[iRow],value);
5397        rowObjectiveWork_[iRow] += value;
5398      }
5399    }
5400  }
5401  // 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};
5402  // 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};
5403  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};
5404  //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};
5405  double extraWeight=10.0;
5406  // Scale back if wanted
5407  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};
5408  if (constantPerturbation<99.0*dualTolerance_) {
5409    perturbation *= 0.1;
5410    extraWeight=0.5;
5411    memcpy(weight,weight2,sizeof(weight2));
5412  }
5413  // adjust weights if all columns long
5414  double factor=1.0;
5415  if (maxLength) {
5416    factor = 3.0/static_cast<double> (minLength);
5417  }
5418  // Make variables with more elements more expensive
5419  const double m1 = 0.5;
5420  double smallestAllowed = CoinMin(1.0e-2*dualTolerance_,maximumFraction);
5421  double largestAllowed = CoinMax(1.0e3*dualTolerance_,maximumFraction*averageCost);
5422  // smaller if in BAB
5423  //if (inCbcOrOther)
5424  //largestAllowed=CoinMin(largestAllowed,1.0e-5);
5425  //smallestAllowed = CoinMin(smallestAllowed,0.1*largestAllowed);
5426#define SAVE_PERT
5427#ifdef SAVE_PERT
5428  if (2*numberColumns_>maximumPerturbationSize_) {
5429    delete [] perturbationArray_;
5430    maximumPerturbationSize_ = 2* numberColumns_;
5431    perturbationArray_ = new double [maximumPerturbationSize_];
5432    for (iColumn=0;iColumn<maximumPerturbationSize_;iColumn++) {
5433      perturbationArray_[iColumn] = randomNumberGenerator_.randomDouble();
5434    }
5435  }
5436#endif
5437  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
5438    if (columnLowerWork_[iColumn]<columnUpperWork_[iColumn]&&getStatus(iColumn)!=basic) {
5439      double value = perturbation;
5440      double currentValue = objectiveWork_[iColumn];
5441      value = CoinMin(value,constantPerturbation+maximumFraction*(fabs(currentValue)+1.0e-1*perturbation+1.0e-8));
5442      //value = CoinMin(value,constantPerturbation;+maximumFraction*fabs(currentValue));
5443      double value2 = constantPerturbation+1.0e-1*smallestNonZero;
5444      if (uniformChange) {
5445        value = maximumFraction;
5446        value2=maximumFraction;
5447      }
5448      if (columnLowerWork_[iColumn]>-largeValue_) {
5449        if (fabs(columnLowerWork_[iColumn])<
5450            fabs(columnUpperWork_[iColumn])) {
5451#ifndef SAVE_PERT
5452          value *= (1.0-m1 +m1*randomNumberGenerator_.randomDouble());
5453          value2 *= (1.0-m1 +m1*randomNumberGenerator_.randomDouble());
5454#else
5455          value *= (1.0-m1 +m1*perturbationArray_[2*iColumn]);
5456          value2 *= (1.0-m1 +m1*perturbationArray_[2*iColumn+1]);
5457#endif
5458        } else {
5459          //value *= -(1.0-m1+m1*randomNumberGenerator_.randomDouble());
5460          //value2 *= -(1.0-m1+m1*randomNumberGenerator_.randomDouble());
5461          value=0.0;
5462        }
5463      } else if (columnUpperWork_[iColumn]<largeValue_) {
5464#ifndef SAVE_PERT
5465        value *= -(1.0-m1+m1*randomNumberGenerator_.randomDouble());
5466        value2 *= -(1.0-m1+m1*randomNumberGenerator_.randomDouble());
5467#else
5468        value *= -(1.0-m1+m1*perturbationArray_[2*iColumn]);
5469        value2 *= -(1.0-m1+m1*perturbationArray_[2*iColumn+1]);
5470#endif
5471      } else {
5472        value=0.0;
5473      }
5474      if (value) {
5475        int length = columnLength[iColumn];
5476        if (length>3) {
5477          length = static_cast<int> (static_cast<double> (length) * factor);
5478          length = CoinMax(3,length);
5479        }
5480        double multiplier;
5481#if 1
5482        if (length<10)
5483          multiplier=weight[length];
5484        else
5485          multiplier = weight[10];
5486#else
5487        if (length<10)
5488          multiplier=weight[length];
5489        else
5490          multiplier = weight[10]+extraWeight*(length-10);
5491        multiplier *= 0.5;
5492#endif
5493        value *= multiplier;
5494        value = CoinMin(value,value2);
5495        if (savePerturbation<50||savePerturbation>60) {
5496          if (fabs(value)<=dualTolerance_)
5497            value=0.0;
5498        } else if (value) {
5499          // get in range
5500          if (fabs(value)<=smallestAllowed) {
5501            value *= 10.0;
5502            while (fabs(value)<=smallestAllowed) 
5503              value *= 10.0;
5504          } else if (fabs(value)>largestAllowed) {
5505            value *= 0.1;
5506            while (fabs(value)>largestAllowed) 
5507              value *= 0.1;
5508          }
5509        }
5510        if (currentValue) {
5511          largest = CoinMax(largest,fabs(value));
5512          if (fabs(value)>fabs(currentValue)*largestPerCent)
5513            largestPerCent=fabs(value/currentValue);
5514        } else {
5515          largestZero=CoinMax(largestZero,fabs(value));
5516        }
5517        // but negative if at ub
5518        if (getStatus(iColumn)==atUpperBound)
5519          value = -value;
5520        if (printOut)
5521          printf("col %d cost %g change %g\n",iColumn,objectiveWork_[iColumn],value);
5522        objectiveWork_[iColumn] += value;
5523      }
5524    }
5525  }
5526  handler_->message(CLP_SIMPLEX_PERTURB,messages_)
5527    <<100.0*maximumFraction<<perturbation<<largest<<100.0*largestPerCent<<largestZero
5528    <<CoinMessageEol;
5529  // and zero changes
5530  //int nTotal = numberRows_+numberColumns_;
5531  //CoinZeroN(cost_+nTotal,nTotal);
5532  // say perturbed
5533  perturbation_=101;
5534  return 0;
5535}
5536/* For strong branching.  On input lower and upper are new bounds
5537   while on output they are change in objective function values
5538   (>1.0e50 infeasible).
5539   Return code is 0 if nothing interesting, -1 if infeasible both
5540   ways and +1 if infeasible one way (check values to see which one(s))
5541   Returns -2 if bad factorization
5542*/
5543int ClpSimplexDual::strongBranching(int numberVariables,const int * variables,
5544                                    double * newLower, double * newUpper,
5545                                    double ** outputSolution,
5546                                    int * outputStatus, int * outputIterations,
5547                                    bool stopOnFirstInfeasible,
5548                                    bool alwaysFinish,
5549                                    int startFinishOptions)
5550{
5551  int i;
5552  int returnCode=0;
5553  double saveObjectiveValue = objectiveValue_;
5554  algorithm_ = -1;
5555
5556  //scaling(false);
5557
5558  // put in standard form (and make row copy)
5559  // create modifiable copies of model rim and do optional scaling
5560  createRim(7+8+16+32,true,startFinishOptions);
5561
5562  // change newLower and newUpper if scaled
5563
5564  // Do initial factorization
5565  // and set certain stuff
5566  // We can either set increasing rows so ...IsBasic gives pivot row
5567  // or we can just increment iBasic one by one
5568  // for now let ...iBasic give pivot row
5569  int useFactorization=false;
5570  if ((startFinishOptions&2)!=0&&(whatsChanged_&(2+512))==2+512)
5571    useFactorization=true; // Keep factorization if possible
5572  // switch off factorization if bad
5573  if (pivotVariable_[0]<0)
5574    useFactorization=false;
5575  if (!useFactorization||factorization_->numberRows()!=numberRows_) {
5576    useFactorization = false;
5577    factorization_->setDefaultValues();
5578
5579    int factorizationStatus = internalFactorize(0);
5580    if (factorizationStatus<0) {
5581      // some error
5582      // we should either debug or ignore
5583#ifndef NDEBUG
5584      printf("***** ClpDual strong branching factorization error - debug\n");
5585#endif
5586      return -2;
5587    } else if (factorizationStatus&&factorizationStatus<=numberRows_) {
5588      handler_->message(CLP_SINGULARITIES,messages_)
5589        <<factorizationStatus
5590        <<CoinMessageEol;
5591    }
5592  }
5593  // save stuff
5594  ClpFactorization saveFactorization(*factorization_);
5595  // save basis and solution
5596  double * saveSolution = new double[numberRows_+numberColumns_];
5597  CoinMemcpyN(solution_,
5598         numberRows_+numberColumns_,saveSolution);
5599  unsigned char * saveStatus = 
5600    new unsigned char [numberRows_+numberColumns_];
5601  CoinMemcpyN(status_,numberColumns_+numberRows_,saveStatus);
5602  // save bounds as createRim makes clean copies
5603  double * saveLower = new double[numberRows_+numberColumns_];
5604  CoinMemcpyN(lower_,
5605         numberRows_+numberColumns_,saveLower);
5606  double * saveUpper = new double[numberRows_+numberColumns_];
5607  CoinMemcpyN(upper_,
5608         numberRows_+numberColumns_,saveUpper);
5609  double * saveObjective = new double[numberRows_+numberColumns_];
5610  CoinMemcpyN(cost_,
5611         numberRows_+numberColumns_,saveObjective);
5612  int * savePivot = new int [numberRows_];
5613  CoinMemcpyN(pivotVariable_, numberRows_,savePivot);
5614  // need to save/restore weights.
5615
5616  int iSolution = 0;
5617  for (i=0;i<numberVariables;i++) {
5618    int iColumn = variables[i];
5619    double objectiveChange;
5620    double saveBound;
5621   
5622    // try down
5623   
5624    saveBound = columnUpper_[iColumn];
5625    // external view - in case really getting optimal
5626    columnUpper_[iColumn] = newUpper[i];
5627    if (scalingFlag_<=0) 
5628      upper_[iColumn] = newUpper[i]*rhsScale_;
5629    else 
5630      upper_[iColumn] = (newUpper[i]/columnScale_[iColumn])*rhsScale_; // scale
5631    // Get fake bounds correctly
5632    double changeCost;
5633    changeBounds(3,NULL,changeCost);
5634    // Start of fast iterations
5635    int status = fastDual(alwaysFinish);
5636    CoinAssert (problemStatus_||objectiveValue_<1.0e50);
5637    // make sure plausible
5638    double obj = CoinMax(objectiveValue_,saveObjectiveValue);
5639    if (status&&problemStatus_!=3) {
5640      // not finished - might be optimal
5641      checkPrimalSolution(rowActivityWork_,columnActivityWork_);
5642      double limit = 0.0;
5643      getDblParam(ClpDualObjectiveLimit, limit);
5644      if (!numberPrimalInfeasibilities_&&obj<limit) { 
5645        problemStatus_=0;
5646      } 
5647      status=problemStatus_;
5648    }
5649    if (status||(problemStatus_==0&&!isDualObjectiveLimitReached())) {
5650      objectiveChange = obj-saveObjectiveValue;
5651    } else {
5652      objectiveChange = 1.0e100;
5653      status=1;
5654    }
5655    if (problemStatus_==3)
5656      status=2;
5657
5658    if (scalingFlag_<=0) {
5659      CoinMemcpyN(solution_,numberColumns_,outputSolution[iSolution]);
5660    } else {
5661      int j;
5662      double * sol = outputSolution[iSolution];
5663      for (j=0;j<numberColumns_;j++) 
5664        sol[j] = solution_[j]*columnScale_[j];
5665    }
5666    outputStatus[iSolution]=status;
5667    outputIterations[iSolution]=numberIterations_;
5668    iSolution++;
5669    // restore
5670    CoinMemcpyN(saveSolution,
5671            numberRows_+numberColumns_,solution_);
5672    CoinMemcpyN(saveStatus,numberColumns_+numberRows_,status_);
5673    CoinMemcpyN(saveLower,
5674            numberRows_+numberColumns_,lower_);
5675    CoinMemcpyN(saveUpper,
5676            numberRows_+numberColumns_,upper_);
5677    CoinMemcpyN(saveObjective,
5678            numberRows_+numberColumns_,cost_);
5679    columnUpper_[iColumn] = saveBound;
5680    CoinMemcpyN(savePivot, numberRows_,pivotVariable_);
5681    delete factorization_;
5682    factorization_ = new ClpFactorization(saveFactorization,numberRows_);
5683
5684    newUpper[i]=objectiveChange;
5685#ifdef CLP_DEBUG
5686    printf("down on %d costs %g\n",iColumn,objectiveChange);
5687#endif
5688         
5689    // try up
5690   
5691    saveBound = columnLower_[iColumn];
5692    // external view - in case really getting optimal
5693    columnLower_[iColumn] = newLower[i];
5694    if (scalingFlag_<=0) 
5695      lower_[iColumn] = newLower[i]*rhsScale_;
5696    else 
5697      lower_[iColumn] = (newLower[i]/columnScale_[iColumn])*rhsScale_; // scale
5698    // Get fake bounds correctly
5699    resetFakeBounds(1);
5700    // Start of fast iterations
5701    status = fastDual(alwaysFinish);
5702    // make sure plausible
5703    obj = CoinMax(objectiveValue_,saveObjectiveValue);
5704    if (status&&problemStatus_!=3) {
5705      // not finished - might be optimal
5706      checkPrimalSolution(rowActivityWork_,columnActivityWork_);
5707      double limit = 0.0;
5708      getDblParam(ClpDualObjectiveLimit, limit);
5709      if (!numberPrimalInfeasibilities_&&obj< limit) { 
5710        problemStatus_=0;
5711      } 
5712      status=problemStatus_;
5713    }
5714    if (status||(problemStatus_==0&&!isDualObjectiveLimitReached())) {
5715      objectiveChange = obj-saveObjectiveValue;
5716    } else {
5717      objectiveChange = 1.0e100;
5718      status=1;
5719    }
5720    if (problemStatus_==3)
5721      status=2;
5722    if (scalingFlag_<=0) {
5723      CoinMemcpyN(solution_,numberColumns_,outputSolution[iSolution]);
5724    } else {
5725      int j;
5726      double * sol = outputSolution[iSolution];
5727      for (j=0;j<numberColumns_;j++) 
5728        sol[j] = solution_[j]*columnScale_[j];
5729    }
5730    outputStatus[iSolution]=status;
5731    outputIterations[iSolution]=numberIterations_;
5732    iSolution++;
5733
5734    // restore
5735    CoinMemcpyN(saveSolution,
5736            numberRows_+numberColumns_,solution_);
5737    CoinMemcpyN(saveStatus,numberColumns_+numberRows_,status_);
5738    CoinMemcpyN(saveLower,
5739            numberRows_+numberColumns_,lower_);
5740    CoinMemcpyN(saveUpper,
5741            numberRows_+numberColumns_,upper_);
5742    CoinMemcpyN(saveObjective,
5743            numberRows_+numberColumns_,cost_);
5744    columnLower_[iColumn] = saveBound;
5745    CoinMemcpyN(savePivot, numberRows_,pivotVariable_);
5746    delete factorization_;
5747    factorization_ = new ClpFactorization(saveFactorization,numberRows_);
5748
5749    newLower[i]=objectiveChange;
5750#ifdef CLP_DEBUG
5751    printf("up on %d costs %g\n",iColumn,objectiveChange);
5752#endif
5753         
5754    /* Possibilities are:
5755       Both sides feasible - store
5756       Neither side feasible - set objective high and exit if desired
5757       One side feasible - change bounds and resolve
5758    */
5759    if (newUpper[i]<1.0e100) {
5760      if(newLower[i]<1.0e100) {
5761        // feasible - no action
5762      } else {
5763        // up feasible, down infeasible
5764        returnCode=1;
5765        if (stopOnFirstInfeasible)
5766          break;
5767      }
5768    } else {
5769      if(newLower[i]<1.0e100) {
5770        // down feasible, up infeasible
5771        returnCode=1;
5772        if (stopOnFirstInfeasible)
5773          break;
5774      } else {
5775        // neither side feasible
5776        returnCode=-1;
5777        break;
5778      }
5779    }
5780  }
5781  delete [] saveSolution;
5782  delete [] saveLower;
5783  delete [] saveUpper;
5784  delete [] saveObjective;
5785  delete [] saveStatus;
5786  delete [] savePivot;
5787  if ((startFinishOptions&1)==0) {
5788    deleteRim(1);
5789    whatsChanged_=0;
5790  } else {
5791    // Original factorization will have been put back by last loop
5792    //delete factorization_;
5793    //factorization_ = new ClpFactorization(saveFactorization);
5794    deleteRim(0);
5795    // mark all as current
5796    whatsChanged_ = 0xffff;
5797  }
5798  objectiveValue_ = saveObjectiveValue;
5799  return returnCode;
5800}
5801// treat no pivot as finished (unless interesting)