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

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

take out message

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