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

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

many changes to try and improve performance

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