source: stable/1.8/Clp/src/ClpSimplexDual.cpp @ 1249

Last change on this file since 1249 was 1249, checked in by ladanyi, 13 years ago

Updated stable/1.8 to match trunk

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