source: branches/devel/Clp/src/ClpSimplexDual.cpp @ 989

Last change on this file since 989 was 989, checked in by forrest, 13 years ago

this may be a mistake

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