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

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

moving branches/devel to trunk

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 174.9 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    // But should not move objectivetoo much ??
3191#define DONT_MOVE_OBJECTIVE
3192#ifdef DONT_MOVE_OBJECTIVE
3193    double moveObjective = fabs(modification*solution_[sequenceIn_]);
3194    double smallMove = CoinMax(fabs(objectiveValue_),1.0e-3);
3195    if (moveObjective>smallMove) {
3196      printf("would move objective by %g - original mod %g sol value %g\n",moveObjective,
3197             modification,solution_[sequenceIn_]);
3198      modification *= smallMove/moveObjective;
3199    }
3200#endif
3201    if ((specialOptions_&(2048+4096))!=0) {
3202      if ((specialOptions_&16384)!=0) {
3203        // in fast dual
3204        if (fabs(modification)<1.0e-7)
3205          modification=0.0;
3206      } else if ((specialOptions_&2048)!=0) {
3207        if (fabs(modification)<1.0e-10)
3208          modification=0.0;
3209      } else {
3210        if (fabs(modification)<1.0e-12)
3211          modification=0.0;
3212      }
3213    }
3214    dualIn_ += modification;
3215    dj_[sequenceIn_]=dualIn_;
3216    cost_[sequenceIn_] += modification;
3217    if (modification)
3218      numberChanged_ ++; // Say changed costs
3219    //int costOffset = numberRows_+numberColumns_;
3220    //cost_[sequenceIn_+costOffset] += modification; // save change
3221    //assert (fabs(modification)<1.0e-6);
3222#ifdef CLP_DEBUG
3223    if ((handler_->logLevel()&32)&&fabs(modification)>1.0e-15)
3224      printf("exact %d new cost %g, change %g\n",sequenceIn_,
3225             cost_[sequenceIn_],modification);
3226#endif
3227#endif
3228   
3229    if (alpha_<0.0) {
3230      // as if from upper bound
3231      directionIn_=-1;
3232      upperIn_=valueIn_;
3233    } else {
3234      // as if from lower bound
3235      directionIn_=1;
3236      lowerIn_=valueIn_;
3237    }
3238  }
3239  //if (thisIncrease)
3240  //dualTolerance_+= 1.0e-6*dblParam_[ClpDualTolerance];
3241
3242  // clear arrays
3243
3244  for (i=0;i<2;i++) {
3245    CoinZeroN(array[i],marker[i][0]);
3246    CoinZeroN(array[i]+marker[i][1],numberColumns_-marker[i][1]);
3247  }
3248  return bestPossible;
3249}
3250#ifdef CLP_ALL_ONE_FILE
3251#undef MAXTRY
3252#endif
3253/* Checks if tentative optimal actually means unbounded
3254   Returns -3 if not, 2 if is unbounded */
3255int 
3256ClpSimplexDual::checkUnbounded(CoinIndexedVector * ray,
3257                                   CoinIndexedVector * spare,
3258                                   double changeCost)
3259{
3260  int status=2; // say unbounded
3261  factorization_->updateColumn(spare,ray);
3262  // get reduced cost
3263  int i;
3264  int number=ray->getNumElements();
3265  int * index = ray->getIndices();
3266  double * array = ray->denseVector();
3267  for (i=0;i<number;i++) {
3268    int iRow=index[i];
3269    int iPivot=pivotVariable_[iRow];
3270    changeCost -= cost(iPivot)*array[iRow];
3271  }
3272  double way;
3273  if (changeCost>0.0) {
3274    //try going down
3275    way=1.0;
3276  } else if (changeCost<0.0) {
3277    //try going up
3278    way=-1.0;
3279  } else {
3280#ifdef CLP_DEBUG
3281    printf("can't decide on up or down\n");
3282#endif
3283    way=0.0;
3284    status=-3;
3285  }
3286  double movement=1.0e10*way; // some largish number
3287  double zeroTolerance = 1.0e-14*dualBound_;
3288  for (i=0;i<number;i++) {
3289    int iRow=index[i];
3290    int iPivot=pivotVariable_[iRow];
3291    double arrayValue = array[iRow];
3292    if (fabs(arrayValue)<zeroTolerance)
3293      arrayValue=0.0;
3294    double newValue=solution(iPivot)+movement*arrayValue;
3295    if (newValue>upper(iPivot)+primalTolerance_||
3296        newValue<lower(iPivot)-primalTolerance_)
3297      status=-3; // not unbounded
3298  }
3299  if (status==2) {
3300    // create ray
3301    delete [] ray_;
3302    ray_ = new double [numberColumns_];
3303    CoinZeroN(ray_,numberColumns_);
3304    for (i=0;i<number;i++) {
3305      int iRow=index[i];
3306      int iPivot=pivotVariable_[iRow];
3307      double arrayValue = array[iRow];
3308      if (iPivot<numberColumns_&&fabs(arrayValue)>=zeroTolerance)
3309        ray_[iPivot] = way* array[iRow];
3310    }
3311  }
3312  ray->clear();
3313  return status;
3314}
3315//static int count_alpha=0;
3316/* Checks if finished.  Updates status */
3317void 
3318ClpSimplexDual::statusOfProblemInDual(int & lastCleaned,int type,
3319                                      double * givenDuals, ClpDataSave & saveData,
3320                                      int ifValuesPass)
3321{
3322  // If lots of iterations then adjust costs if large ones
3323  if (numberIterations_>4*(numberRows_+numberColumns_)&&objectiveScale_==1.0) {
3324    double largest=0.0;
3325    for (int i=0;i<numberRows_;i++) {
3326      int iColumn = pivotVariable_[i];
3327      largest = CoinMax(largest,fabs(cost_[iColumn]));
3328    }
3329    if (largest>1.0e6) {
3330      objectiveScale_ = 1.0e6/largest;
3331      for (int i=0;i<numberRows_+numberColumns_;i++)
3332        cost_[i] *= objectiveScale_;
3333    }
3334  }
3335  bool normalType=true;
3336  int numberPivots = factorization_->pivots();
3337  double realDualInfeasibilities=0.0;
3338  if (type==2) {
3339    if (alphaAccuracy_!=-1.0)
3340      alphaAccuracy_=-2.0;
3341    // trouble - restore solution
3342    CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
3343    CoinMemcpyN(savedSolution_+numberColumns_ ,
3344           numberRows_,rowActivityWork_);
3345    CoinMemcpyN(savedSolution_ ,
3346           numberColumns_,columnActivityWork_);
3347    // restore extra stuff
3348    int dummy;
3349    matrix_->generalExpanded(this,6,dummy);
3350    forceFactorization_=1; // a bit drastic but ..
3351    changeMade_++; // say something changed
3352    // get correct bounds on all variables
3353    resetFakeBounds();
3354  }
3355  int tentativeStatus = problemStatus_;
3356  double changeCost;
3357  bool unflagVariables = true;
3358  bool weightsSaved=false;
3359  if (alphaAccuracy_<0.0||!numberPivots||alphaAccuracy_>1.0e4||factorization_->pivots()>20) {
3360    if (problemStatus_>-3||numberPivots) {
3361      // factorize
3362      // later on we will need to recover from singularities
3363      // also we could skip if first time
3364      // save dual weights
3365      dualRowPivot_->saveWeights(this,1);
3366      weightsSaved=true;
3367      if (type) {
3368        // is factorization okay?
3369        if (internalFactorize(1)) {
3370          // no - restore previous basis
3371          unflagVariables = false;
3372          assert (type==1);
3373          changeMade_++; // say something changed
3374          // Keep any flagged variables
3375          int i;
3376          for (i=0;i<numberRows_+numberColumns_;i++) {
3377            if (flagged(i))
3378              saveStatus_[i] |= 64; //say flagged
3379          }
3380          CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
3381          CoinMemcpyN(savedSolution_+numberColumns_ ,
3382                      numberRows_,rowActivityWork_);
3383          CoinMemcpyN(savedSolution_ ,
3384                      numberColumns_,columnActivityWork_);
3385          // restore extra stuff
3386          int dummy;
3387          matrix_->generalExpanded(this,6,dummy);
3388          // get correct bounds on all variables
3389          resetFakeBounds();
3390          // need to reject something
3391          char x = isColumn(sequenceOut_) ? 'C' :'R';
3392          handler_->message(CLP_SIMPLEX_FLAG,messages_)
3393            <<x<<sequenceWithin(sequenceOut_)
3394            <<CoinMessageEol;
3395          setFlagged(sequenceOut_);
3396          progress_->clearBadTimes();
3397         
3398          // Go to safe
3399          factorization_->pivotTolerance(0.99);
3400          forceFactorization_=1; // a bit drastic but ..
3401          type = 2;
3402          //assert (internalFactorize(1)==0);
3403          if (internalFactorize(1)) {
3404            CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
3405            CoinMemcpyN(savedSolution_+numberColumns_ ,
3406                        numberRows_,rowActivityWork_);
3407            CoinMemcpyN(savedSolution_ ,
3408                        numberColumns_,columnActivityWork_);
3409            // restore extra stuff
3410            int dummy;
3411            matrix_->generalExpanded(this,6,dummy);
3412            // debug
3413            int returnCode = internalFactorize(1);
3414            while (returnCode) {
3415              // ouch
3416              // switch off dense
3417              int saveDense = factorization_->denseThreshold();
3418              factorization_->setDenseThreshold(0);
3419              // Go to safe
3420              factorization_->pivotTolerance(0.99);
3421              // make sure will do safe factorization
3422              pivotVariable_[0]=-1;
3423              returnCode=internalFactorize(2);
3424              factorization_->setDenseThreshold(saveDense);
3425            }
3426            // get correct bounds on all variables
3427            resetFakeBounds();
3428          }
3429        }
3430      }
3431      if (problemStatus_!=-4||numberPivots>10)
3432        problemStatus_=-3;
3433    }
3434  } else {
3435    //printf("testing with accuracy of %g and status of %d\n",alphaAccuracy_,problemStatus_);
3436    //count_alpha++;
3437    //if ((count_alpha%5000)==0)
3438    //printf("count alpha %d\n",count_alpha);
3439  }
3440  // at this stage status is -3 or -4 if looks infeasible
3441  // get primal and dual solutions
3442  gutsOfSolution(givenDuals,NULL);
3443  // If bad accuracy treat as singular
3444  if ((largestPrimalError_>1.0e15||largestDualError_>1.0e15)&&numberIterations_) {
3445    // restore previous basis
3446    unflagVariables = false;
3447    changeMade_++; // say something changed
3448    // Keep any flagged variables
3449    int i;
3450    for (i=0;i<numberRows_+numberColumns_;i++) {
3451      if (flagged(i))
3452        saveStatus_[i] |= 64; //say flagged
3453    }
3454    CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
3455    CoinMemcpyN(savedSolution_+numberColumns_ ,
3456           numberRows_,rowActivityWork_);
3457    CoinMemcpyN(savedSolution_ ,
3458           numberColumns_,columnActivityWork_);
3459    // restore extra stuff
3460    int dummy;
3461    matrix_->generalExpanded(this,6,dummy);
3462    // get correct bounds on all variables
3463    resetFakeBounds();
3464    // need to reject something
3465    char x = isColumn(sequenceOut_) ? 'C' :'R';
3466    handler_->message(CLP_SIMPLEX_FLAG,messages_)
3467      <<x<<sequenceWithin(sequenceOut_)
3468      <<CoinMessageEol;
3469    setFlagged(sequenceOut_);
3470    progress_->clearBadTimes();
3471   
3472    // Go to safer
3473    double newTolerance = CoinMin(1.1*factorization_->pivotTolerance(),0.99);
3474    factorization_->pivotTolerance(newTolerance);
3475    forceFactorization_=1; // a bit drastic but ..
3476    if (alphaAccuracy_!=-1.0)
3477      alphaAccuracy_=-2.0;
3478    type = 2;
3479    //assert (internalFactorize(1)==0);
3480    if (internalFactorize(1)) {
3481      CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
3482      CoinMemcpyN(savedSolution_+numberColumns_ ,
3483              numberRows_,rowActivityWork_);
3484      CoinMemcpyN(savedSolution_ ,
3485           numberColumns_,columnActivityWork_);
3486      // restore extra stuff
3487      int dummy;
3488      matrix_->generalExpanded(this,6,dummy);
3489      // debug
3490      int returnCode = internalFactorize(1);
3491      while (returnCode) {
3492        // ouch
3493        // switch off dense
3494        int saveDense = factorization_->denseThreshold();
3495        factorization_->setDenseThreshold(0);
3496        // Go to safe
3497        factorization_->pivotTolerance(0.99);
3498        // make sure will do safe factorization
3499        pivotVariable_[0]=-1;
3500        returnCode=internalFactorize(2);
3501        factorization_->setDenseThreshold(saveDense);
3502      }
3503      // get correct bounds on all variables
3504      resetFakeBounds();
3505    }
3506    // get primal and dual solutions
3507    gutsOfSolution(givenDuals,NULL);
3508  } else if (largestPrimalError_<1.0e-7&&largestDualError_<1.0e-7) {
3509    // Can reduce tolerance
3510    double newTolerance = CoinMax(0.99*factorization_->pivotTolerance(),saveData.pivotTolerance_);
3511    factorization_->pivotTolerance(newTolerance);
3512  } 
3513  // Double check infeasibility if no action
3514  if (progress_->lastIterationNumber(0)==numberIterations_) {
3515    if (dualRowPivot_->looksOptimal()) {
3516      numberPrimalInfeasibilities_ = 0;
3517      sumPrimalInfeasibilities_ = 0.0;
3518    }
3519#if 1
3520  } else {
3521    double thisObj = objectiveValue_;
3522    double lastObj = progress_->lastObjective(0);
3523    if (lastObj>thisObj+1.0e-3*CoinMax(fabs(thisObj),fabs(lastObj))+1.0
3524        &&!ifValuesPass&&firstFree_<0) {
3525      int maxFactor = factorization_->maximumPivots();
3526      if (maxFactor>10&&numberPivots>1) {
3527        //printf("lastobj %g thisobj %g\n",lastObj,thisObj);
3528        //if (forceFactorization_<0)
3529        //forceFactorization_= maxFactor;
3530        //forceFactorization_ = CoinMax(1,(forceFactorization_>>1));
3531        forceFactorization_=1;
3532        //printf("Reducing factorization frequency - bad backwards\n");
3533        unflagVariables = false;
3534        changeMade_++; // say something changed
3535        CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
3536        CoinMemcpyN(savedSolution_+numberColumns_ ,
3537                numberRows_,rowActivityWork_);
3538        CoinMemcpyN(savedSolution_ ,
3539                numberColumns_,columnActivityWork_);
3540        // restore extra stuff
3541        int dummy;
3542        matrix_->generalExpanded(this,6,dummy);
3543        // get correct bounds on all variables
3544        resetFakeBounds();
3545        if(factorization_->pivotTolerance()<0.2)
3546          factorization_->pivotTolerance(0.2);
3547        if (alphaAccuracy_!=-1.0)
3548          alphaAccuracy_=-2.0;
3549        if (internalFactorize(1)) {
3550          CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
3551          CoinMemcpyN(savedSolution_+numberColumns_ ,
3552                  numberRows_,rowActivityWork_);
3553          CoinMemcpyN(savedSolution_ ,
3554                  numberColumns_,columnActivityWork_);
3555          // restore extra stuff
3556          int dummy;
3557          matrix_->generalExpanded(this,6,dummy);
3558          // debug
3559          int returnCode = internalFactorize(1);
3560          while (returnCode) {
3561            // ouch
3562            // switch off dense
3563            int saveDense = factorization_->denseThreshold();
3564            factorization_->setDenseThreshold(0);
3565            // Go to safe
3566            factorization_->pivotTolerance(0.99);
3567            // make sure will do safe factorization
3568            pivotVariable_[0]=-1;
3569            returnCode=internalFactorize(2);
3570            factorization_->setDenseThreshold(saveDense);
3571          }
3572          resetFakeBounds();
3573        }
3574        type = 2; // so will restore weights
3575        // get primal and dual solutions
3576        gutsOfSolution(givenDuals,NULL);
3577      } 
3578    }
3579#endif
3580  }
3581  // Up tolerance if looks a bit odd
3582  if (numberIterations_>CoinMax(1000,numberRows_>>4)&&(specialOptions_&64)!=0) {
3583    if (sumPrimalInfeasibilities_&&sumPrimalInfeasibilities_<1.0e5) {
3584      int backIteration = progress_->lastIterationNumber(CLP_PROGRESS-1);
3585      if (backIteration>0&&numberIterations_-backIteration<9*CLP_PROGRESS) {
3586        if (factorization_->pivotTolerance()<0.9) {
3587          // up tolerance
3588          factorization_->pivotTolerance(CoinMin(factorization_->pivotTolerance()*1.05+0.02,0.91));
3589          //printf("tol now %g\n",factorization_->pivotTolerance());
3590          progress_->clearIterationNumbers();
3591        }
3592      }
3593    }
3594  }
3595  // Check if looping
3596  int loop;
3597  if (!givenDuals&&type!=2) 
3598    loop = progress_->looping();
3599  else
3600    loop=-1;
3601  int situationChanged=0;
3602  if (loop>=0) {
3603    problemStatus_ = loop; //exit if in loop
3604    if (!problemStatus_) {
3605      // declaring victory
3606      numberPrimalInfeasibilities_ = 0;
3607      sumPrimalInfeasibilities_ = 0.0;
3608    } else {
3609      problemStatus_ = 10; // instead - try other algorithm
3610    }
3611    return;
3612  } else if (loop<-1) {
3613    // something may have changed
3614    gutsOfSolution(NULL,NULL);
3615    situationChanged=1;
3616  }
3617  // really for free variables in
3618  if((progressFlag_&2)!=0) {
3619    situationChanged=2;
3620  }
3621  progressFlag_ = 0; //reset progress flag
3622  if (handler_->detail(CLP_SIMPLEX_STATUS,messages_)<100) {
3623    handler_->message(CLP_SIMPLEX_STATUS,messages_)
3624      <<numberIterations_<<objectiveValue();
3625    handler_->printing(sumPrimalInfeasibilities_>0.0)
3626      <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
3627    handler_->printing(sumDualInfeasibilities_>0.0)
3628      <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
3629    handler_->printing(numberDualInfeasibilitiesWithoutFree_
3630                       <numberDualInfeasibilities_)
3631                         <<numberDualInfeasibilitiesWithoutFree_;
3632    handler_->message()<<CoinMessageEol;
3633  }
3634  realDualInfeasibilities=sumDualInfeasibilities_;
3635  double saveTolerance =dualTolerance_;
3636  // If we need to carry on cleaning variables
3637  if (!numberPrimalInfeasibilities_&&(specialOptions_&1024)!=0&&CLEAN_FIXED) {
3638    for (int iRow=0;iRow<numberRows_;iRow++) {
3639      int iPivot=pivotVariable_[iRow];
3640      if (!flagged(iPivot)&&pivoted(iPivot)) {
3641        // carry on
3642        numberPrimalInfeasibilities_=-1;
3643        sumOfRelaxedPrimalInfeasibilities_ = 1.0;
3644        sumPrimalInfeasibilities_ = 1.0;
3645        break;
3646      }
3647    }
3648  }
3649  /* If we are primal feasible and any dual infeasibilities are on
3650     free variables then it is better to go to primal */
3651  if (!numberPrimalInfeasibilities_&&!numberDualInfeasibilitiesWithoutFree_&&
3652      numberDualInfeasibilities_)
3653    problemStatus_=10;
3654  // dual bound coming in
3655  //double saveDualBound = dualBound_;
3656  while (problemStatus_<=-3) {
3657    int cleanDuals=0;
3658    if (situationChanged!=0)
3659      cleanDuals=1;
3660    int numberChangedBounds=0;
3661    int doOriginalTolerance=0;
3662    if ( lastCleaned==numberIterations_)
3663      doOriginalTolerance=1;
3664    // check optimal
3665    // give code benefit of doubt
3666    if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
3667        sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
3668      // say optimal (with these bounds etc)
3669      numberDualInfeasibilities_ = 0;
3670      sumDualInfeasibilities_ = 0.0;
3671      numberPrimalInfeasibilities_ = 0;
3672      sumPrimalInfeasibilities_ = 0.0;
3673    }
3674    //if (dualFeasible()||problemStatus_==-4||(primalFeasible()&&!numberDualInfeasibilitiesWithoutFree_)) {
3675    if (dualFeasible()||problemStatus_==-4) {
3676      progress_->modifyObjective(objectiveValue_
3677                               -sumDualInfeasibilities_*dualBound_);
3678      if (primalFeasible()&&!givenDuals) {
3679        normalType=false;
3680        // may be optimal - or may be bounds are wrong
3681        handler_->message(CLP_DUAL_BOUNDS,messages_)
3682          <<dualBound_
3683          <<CoinMessageEol;
3684        // save solution in case unbounded
3685        double * saveColumnSolution = NULL;
3686        double * saveRowSolution = NULL;
3687        bool inCbc = (specialOptions_&0x01000000)!=0;
3688        if (!inCbc) {
3689          saveColumnSolution = CoinCopyOfArray(columnActivityWork_,numberColumns_);
3690          saveRowSolution = CoinCopyOfArray(rowActivityWork_,numberRows_);
3691        }
3692        numberChangedBounds=changeBounds(false,rowArray_[3],changeCost);
3693        if (numberChangedBounds<=0&&!numberDualInfeasibilities_) {
3694          //looks optimal - do we need to reset tolerance
3695          if (perturbation_==101) {
3696            perturbation_=102; // stop any perturbations
3697            cleanDuals=1;
3698            // make sure fake bounds are back
3699            //computeObjectiveValue();
3700            changeBounds(true,NULL,changeCost);
3701            //computeObjectiveValue();
3702            createRim(4);
3703            // make sure duals are current
3704            computeDuals(givenDuals);
3705            checkDualSolution();
3706            //if (numberDualInfeasibilities_)
3707              numberChanged_=1; // force something to happen
3708            //else
3709            //computeObjectiveValue();
3710          }
3711          if (lastCleaned<numberIterations_&&numberTimesOptimal_<4&&
3712              (numberChanged_||(specialOptions_&4096)==0)) {
3713            doOriginalTolerance=2;
3714            numberTimesOptimal_++;
3715            changeMade_++; // say something changed
3716            if (numberTimesOptimal_==1) {
3717              dualTolerance_ = dblParam_[ClpDualTolerance];
3718            } else {
3719              if (numberTimesOptimal_==2) {
3720                // better to have small tolerance even if slower
3721                factorization_->zeroTolerance(1.0e-15);
3722              }
3723              dualTolerance_ = dblParam_[ClpDualTolerance];
3724              dualTolerance_ *= pow(2.0,numberTimesOptimal_-1);
3725            }
3726            cleanDuals=2; // If nothing changed optimal else primal
3727          } else {
3728            problemStatus_=0; // optimal
3729            if (lastCleaned<numberIterations_&&numberChanged_) {
3730              handler_->message(CLP_SIMPLEX_GIVINGUP,messages_)
3731                <<CoinMessageEol;
3732            }
3733          }
3734        } else {
3735          cleanDuals=1;
3736          if (doOriginalTolerance==1) {
3737            // check unbounded
3738            // find a variable with bad dj
3739            int iSequence;
3740            int iChosen=-1;
3741            if (!inCbc) {
3742              double largest = 100.0*primalTolerance_;
3743              for (iSequence=0;iSequence<numberRows_+numberColumns_;
3744                   iSequence++) {
3745                double djValue = dj_[iSequence];
3746                double originalLo = originalLower(iSequence);
3747                double originalUp = originalUpper(iSequence);
3748                if (fabs(djValue)>fabs(largest)) {
3749                  if (getStatus(iSequence)!=basic) {
3750                    if (djValue>0&&originalLo<-1.0e20) {
3751                      if (djValue>fabs(largest)) {
3752                        largest=djValue;
3753                        iChosen=iSequence;
3754                      }
3755                    } else if (djValue<0&&originalUp>1.0e20) {
3756                      if (-djValue>fabs(largest)) {
3757                        largest=djValue;
3758                        iChosen=iSequence;
3759                      }
3760                    } 
3761                  }
3762                }
3763              }
3764            }
3765            if (iChosen>=0) {
3766              int iSave=sequenceIn_;
3767              sequenceIn_=iChosen;
3768              unpack(rowArray_[1]);
3769              sequenceIn_ = iSave;
3770              // if dual infeasibilities then must be free vector so add in dual
3771              if (numberDualInfeasibilities_) {
3772                if (fabs(changeCost)>1.0e-5)
3773                  printf("Odd free/unbounded combo\n");
3774                changeCost += cost_[iChosen];
3775              }
3776              problemStatus_ = checkUnbounded(rowArray_[1],rowArray_[0],
3777                                              changeCost);
3778              rowArray_[1]->clear();
3779            } else {
3780              problemStatus_=-3;
3781            }
3782            if (problemStatus_==2&&perturbation_==101) {
3783              perturbation_=102; // stop any perturbations
3784              cleanDuals=1;
3785              createRim(4);
3786              problemStatus_=-1;
3787            }
3788            if (problemStatus_==2) {
3789              // it is unbounded - restore solution
3790              // but first add in changes to non-basic
3791              int iColumn;
3792              double * original = columnArray_[0]->denseVector();
3793              for (iColumn=0;iColumn<numberColumns_;iColumn++) {
3794                if(getColumnStatus(iColumn)!= basic)
3795                  ray_[iColumn] += 
3796                    saveColumnSolution[iColumn]-original[iColumn];
3797                columnActivityWork_[iColumn] = original[iColumn];
3798              }
3799              CoinMemcpyN(saveRowSolution,numberRows_,
3800                                rowActivityWork_);
3801            }
3802          } else {
3803            doOriginalTolerance=2;
3804            rowArray_[0]->clear();
3805          }
3806        }
3807        delete [] saveColumnSolution;
3808        delete [] saveRowSolution;
3809      } 
3810      if (problemStatus_==-4||problemStatus_==-5) {
3811        // may be infeasible - or may be bounds are wrong
3812        numberChangedBounds=changeBounds(false,NULL,changeCost);
3813        /* Should this be here as makes no difference to being feasible.
3814           But seems to make a difference to run times. */
3815        if (perturbation_==101&&0) {
3816          perturbation_=102; // stop any perturbations
3817          cleanDuals=1;
3818          numberChangedBounds=1;
3819          // make sure fake bounds are back
3820          changeBounds(true,NULL,changeCost);
3821          createRim(4);
3822        }
3823        if (numberChangedBounds<=0||dualBound_>1.0e20||
3824            (largestPrimalError_>1.0&&dualBound_>1.0e17)) {
3825          problemStatus_=1; // infeasible
3826          if (perturbation_==101) {
3827            perturbation_=102; // stop any perturbations
3828            //cleanDuals=1;
3829            //numberChangedBounds=1;
3830            //createRim(4);
3831          }
3832        } else {
3833          normalType=false;
3834          problemStatus_=-1; //iterate
3835          cleanDuals=1;
3836          if (numberChangedBounds<=0)
3837            doOriginalTolerance=2;
3838          // and delete ray which has been created
3839          delete [] ray_;
3840          ray_ = NULL;
3841        }
3842
3843      }
3844    } else {
3845      cleanDuals=1;
3846    }
3847    if (problemStatus_<0) {
3848      if (doOriginalTolerance==2) {
3849        // put back original tolerance
3850        lastCleaned=numberIterations_;
3851        numberChanged_ =0; // Number of variables with changed costs
3852        handler_->message(CLP_DUAL_ORIGINAL,messages_)
3853          <<CoinMessageEol;
3854        perturbation_=102; // stop any perturbations
3855#if 0
3856        double * xcost = new double[numberRows_+numberColumns_];
3857        double * xlower = new double[numberRows_+numberColumns_];
3858        double * xupper = new double[numberRows_+numberColumns_];
3859        double * xdj = new double[numberRows_+numberColumns_];
3860        double * xsolution = new double[numberRows_+numberColumns_];
3861        memcpy(xcost,cost_,(numberRows_+numberColumns_)*sizeof(double));
3862        memcpy(xlower,lower_,(numberRows_+numberColumns_)*sizeof(double));
3863        memcpy(xupper,upper_,(numberRows_+numberColumns_)*sizeof(double));
3864        memcpy(xdj,dj_,(numberRows_+numberColumns_)*sizeof(double));
3865        memcpy(xsolution,solution_,(numberRows_+numberColumns_)*sizeof(double));
3866#endif
3867        createRim(4);
3868        // make sure duals are current
3869        computeDuals(givenDuals);
3870        checkDualSolution();
3871#if 0
3872        int i;
3873        for (i=0;i<numberRows_+numberColumns_;i++) {
3874          if (cost_[i]!=xcost[i])
3875            printf("** %d old cost %g new %g sol %g\n",
3876                   i,xcost[i],cost_[i],solution_[i]);
3877          if (lower_[i]!=xlower[i])
3878            printf("** %d old lower %g new %g sol %g\n",
3879                   i,xlower[i],lower_[i],solution_[i]);
3880          if (upper_[i]!=xupper[i])
3881            printf("** %d old upper %g new %g sol %g\n",
3882                   i,xupper[i],upper_[i],solution_[i]);
3883          if (dj_[i]!=xdj[i])
3884            printf("** %d old dj %g new %g sol %g\n",
3885                   i,xdj[i],dj_[i],solution_[i]);
3886          if (solution_[i]!=xsolution[i])
3887            printf("** %d old solution %g new %g sol %g\n",
3888                   i,xsolution[i],solution_[i],solution_[i]);
3889        }
3890        //delete [] xcost;
3891        //delete [] xupper;
3892        //delete [] xlower;
3893        //delete [] xdj;
3894        //delete [] xsolution;
3895#endif
3896        // put back bounds as they were if was optimal
3897        if (doOriginalTolerance==2&&cleanDuals!=2) {
3898          changeMade_++; // say something changed
3899          /* We may have already changed some bounds in this function
3900             so save numberFake_ and add in.
3901
3902             Worst that can happen is that we waste a bit of time  - but it must be finite.
3903          */
3904          int saveNumberFake = numberFake_;
3905          changeBounds(true,NULL,changeCost);
3906          numberFake_ += saveNumberFake;
3907          cleanDuals=2;
3908          //cleanDuals=1;
3909        }
3910#if 0
3911        //int i;
3912        for (i=0;i<numberRows_+numberColumns_;i++) {
3913          if (cost_[i]!=xcost[i])
3914            printf("** %d old cost %g new %g sol %g\n",
3915                   i,xcost[i],cost_[i],solution_[i]);
3916          if (lower_[i]!=xlower[i])
3917            printf("** %d old lower %g new %g sol %g\n",
3918                   i,xlower[i],lower_[i],solution_[i]);
3919          if (upper_[i]!=xupper[i])
3920            printf("** %d old upper %g new %g sol %g\n",
3921                   i,xupper[i],upper_[i],solution_[i]);
3922          if (dj_[i]!=xdj[i])
3923            printf("** %d old dj %g new %g sol %g\n",
3924                   i,xdj[i],dj_[i],solution_[i]);
3925          if (solution_[i]!=xsolution[i])
3926            printf("** %d old solution %g new %g sol %g\n",
3927                   i,xsolution[i],solution_[i],solution_[i]);
3928        }
3929        delete [] xcost;
3930        delete [] xupper;
3931        delete [] xlower;
3932        delete [] xdj;
3933        delete [] xsolution;
3934#endif
3935      }
3936      if (cleanDuals==1||(cleanDuals==2&&!numberDualInfeasibilities_)) {
3937        // make sure dual feasible
3938        // look at all rows and columns
3939        rowArray_[0]->clear();
3940        columnArray_[0]->clear();
3941        double objectiveChange=0.0;
3942#if 0
3943        double * xcost = new double[numberRows_+numberColumns_];
3944        double * xlower = new double[numberRows_+numberColumns_];
3945        double * xupper = new double[numberRows_+numberColumns_];
3946        double * xdj = new double[numberRows_+numberColumns_];
3947        double * xsolution = new double[numberRows_+numberColumns_];
3948        memcpy(xcost,cost_,(numberRows_+numberColumns_)*sizeof(double));
3949        memcpy(xlower,lower_,(numberRows_+numberColumns_)*sizeof(double));
3950        memcpy(xupper,upper_,(numberRows_+numberColumns_)*sizeof(double));
3951        memcpy(xdj,dj_,(numberRows_+numberColumns_)*sizeof(double));
3952        memcpy(xsolution,solution_,(numberRows_+numberColumns_)*sizeof(double));
3953#endif
3954        if (givenDuals)
3955          dualTolerance_=1.0e50;
3956        updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[1],
3957          0.0,objectiveChange,true);
3958        dualTolerance_=saveTolerance;
3959#if 0
3960        int i;
3961        for (i=0;i<numberRows_+numberColumns_;i++) {
3962          if (cost_[i]!=xcost[i])
3963            printf("** %d old cost %g new %g sol %g\n",
3964                   i,xcost[i],cost_[i],solution_[i]);
3965          if (lower_[i]!=xlower[i])
3966            printf("** %d old lower %g new %g sol %g\n",
3967                   i,xlower[i],lower_[i],solution_[i]);
3968          if (upper_[i]!=xupper[i])
3969            printf("** %d old upper %g new %g sol %g\n",
3970                   i,xupper[i],upper_[i],solution_[i]);
3971          if (dj_[i]!=xdj[i])
3972            printf("** %d old dj %g new %g sol %g\n",
3973                   i,xdj[i],dj_[i],solution_[i]);
3974          if (solution_[i]!=xsolution[i])
3975            printf("** %d old solution %g new %g sol %g\n",
3976                   i,xsolution[i],solution_[i],solution_[i]);
3977        }
3978        delete [] xcost;
3979        delete [] xupper;
3980        delete [] xlower;
3981        delete [] xdj;
3982        delete [] xsolution;
3983#endif
3984        // for now - recompute all
3985        gutsOfSolution(NULL,NULL);
3986        if (givenDuals)
3987          dualTolerance_=1.0e50;
3988        updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[1],
3989          0.0,objectiveChange,true);
3990        dualTolerance_=saveTolerance;
3991        //assert(numberDualInfeasibilitiesWithoutFree_==0);
3992        if (numberDualInfeasibilities_) {
3993          if (numberPrimalInfeasibilities_||numberPivots)
3994            problemStatus_=-1; // carry on as normal
3995          else
3996            problemStatus_=10; // try primal
3997        } else if (situationChanged==2) {
3998          problemStatus_=-1; // carry on as normal
3999        }
4000        situationChanged=0;
4001      } else {
4002        // iterate
4003        if (cleanDuals!=2) 
4004          problemStatus_=-1;
4005        else 
4006          problemStatus_ = 10; // try primal
4007      }
4008    }
4009  }
4010  if (type==0||type==1) {
4011    if (!type) {
4012      // create save arrays
4013      delete [] saveStatus_;
4014      delete [] savedSolution_;
4015      saveStatus_ = new unsigned char [numberRows_+numberColumns_];
4016      savedSolution_ = new double [numberRows_+numberColumns_];
4017    }
4018    // save arrays
4019    CoinMemcpyN(status_,numberColumns_+numberRows_,saveStatus_);
4020    CoinMemcpyN(rowActivityWork_,
4021           numberRows_,savedSolution_+numberColumns_);
4022    CoinMemcpyN(columnActivityWork_,numberColumns_,savedSolution_);
4023    // save extra stuff
4024    int dummy;
4025    matrix_->generalExpanded(this,5,dummy);
4026  }
4027  if (weightsSaved) {
4028    // restore weights (if saved) - also recompute infeasibility list
4029    if (tentativeStatus>-3) 
4030      dualRowPivot_->saveWeights(this,(type <2) ? 2 : 4);
4031    else
4032      dualRowPivot_->saveWeights(this,3);
4033  }
4034  // unflag all variables (we may want to wait a bit?)
4035  if ((tentativeStatus!=-2&&tentativeStatus!=-1)&&unflagVariables) {
4036    int iRow;
4037    int numberFlagged=0;
4038    for (iRow=0;iRow<numberRows_;iRow++) {
4039      int iPivot=pivotVariable_[iRow];
4040      if (flagged(iPivot)) {
4041        numberFlagged++;
4042        clearFlagged(iPivot);
4043      }
4044    }
4045#if 0
4046    if (numberFlagged) {
4047      printf("unflagging %d variables - status %d ninf %d nopt %d\n",numberFlagged,tentativeStatus,
4048             numberPrimalInfeasibilities_,
4049             numberTimesOptimal_);
4050    }
4051#endif
4052    unflagVariables = numberFlagged>0;
4053    if (numberFlagged&&!numberPivots) {
4054      /* looks like trouble as we have not done any iterations.
4055         Try changing pivot tolerance then give it a few goes and give up */
4056      if (factorization_->pivotTolerance()<0.9) {
4057        factorization_->pivotTolerance(0.99);
4058      } else if (numberTimesOptimal_<10) {
4059        numberTimesOptimal_++;
4060      } else {
4061        unflagVariables=false;
4062        changeMade_=false;
4063        secondaryStatus_ = 1; // and say probably infeasible
4064      }
4065    }
4066  }
4067  // see if cutoff reached
4068  double limit = 0.0;
4069  getDblParam(ClpDualObjectiveLimit, limit);
4070#if 0
4071  if(fabs(limit)<1.0e30&&objectiveValue()*optimizationDirection_>
4072     limit+100.0) {
4073    printf("lim %g obj %g %g - wo perturb %g sum dual %g\n",
4074           limit,objectiveValue_,objectiveValue(),computeInternalObjectiveValue(),sumDualInfeasibilities_);
4075  }
4076#endif
4077  if(fabs(limit)<1.0e30&&objectiveValue()*optimizationDirection_>
4078     limit&&!numberAtFakeBound()) {
4079    bool looksInfeasible = !numberDualInfeasibilities_;
4080    if (objectiveValue()*optimizationDirection_>limit+fabs(0.1*limit)+1.0e2*sumDualInfeasibilities_+1.0e4&&
4081        sumDualInfeasibilities_<largestDualError_&&numberIterations_>0.5*numberRows_+1000)
4082      looksInfeasible=true;
4083    if (looksInfeasible) {
4084      // Even if not perturbed internal costs may have changed
4085      if (true||perturbation_==101) {
4086        // be careful
4087        if (numberIterations_) {
4088          if(computeInternalObjectiveValue()>limit) {
4089            problemStatus_=1;
4090            secondaryStatus_ = 1; // and say was on cutoff
4091          }
4092        }
4093      } else {
4094        problemStatus_=1;
4095        secondaryStatus_ = 1; // and say was on cutoff
4096      }
4097    }
4098  }
4099  // If we are in trouble and in branch and bound give up
4100  if ((specialOptions_&1024)!=0) {
4101    int looksBad=0;
4102    if (largestPrimalError_*largestDualError_>1.0e2) {
4103      looksBad=1;
4104    } else if (largestPrimalError_>1.0e-2
4105        &&objectiveValue_>CoinMin(1.0e15,1.0e3*limit)) {
4106      looksBad=2;
4107    }
4108    if (looksBad) {
4109      if (factorization_->pivotTolerance()<0.9) {
4110        // up tolerance
4111        factorization_->pivotTolerance(CoinMin(factorization_->pivotTolerance()*1.05+0.02,0.91));
4112      } else if (numberIterations_>10000) {
4113        if (handler_->logLevel()>2)
4114          printf("bad dual - saying infeasible %d\n",looksBad);
4115        problemStatus_=1;
4116        secondaryStatus_ = 1; // and say was on cutoff
4117      } else if (largestPrimalError_>1.0e5) {
4118        {
4119          int iBigB=-1;
4120          double bigB=0.0;
4121          int iBigN=-1;
4122          double bigN=0.0;
4123          for (int i=0;i<numberRows_+numberColumns_;i++) {
4124            double value = fabs(solution_[i]);
4125            if (getStatus(i)==basic) {
4126              if (value>bigB) {
4127                bigB= value;
4128                iBigB=i;
4129              }
4130            } else {
4131              if (value>bigN) {
4132                bigN= value;
4133                iBigN=i;
4134              }
4135            }
4136          }
4137          if (bigB>1.0e8||bigN>1.0e8) {
4138            if (handler_->logLevel()>0)
4139              printf("it %d - basic %d %g, nonbasic %d %g\n",
4140                     numberIterations_,iBigB,bigB,iBigN,bigN);
4141          }
4142        }
4143        if (handler_->logLevel()>2)
4144          printf("bad dual - going to primal %d %g\n",looksBad,largestPrimalError_);
4145        allSlackBasis(true);
4146        problemStatus_=10;
4147      }
4148    }
4149  }
4150  if (problemStatus_<0&&!changeMade_) {
4151    problemStatus_=4; // unknown
4152  }
4153  lastGoodIteration_ = numberIterations_;
4154  if (problemStatus_<0) {
4155    sumDualInfeasibilities_=realDualInfeasibilities; // back to say be careful
4156    if (sumDualInfeasibilities_)
4157      numberDualInfeasibilities_=1;
4158  }
4159#if 1
4160  double thisObj = progress_->lastObjective(0);
4161  double lastObj = progress_->lastObjective(1);
4162  if (lastObj>thisObj+1.0e-4*CoinMax(fabs(thisObj),fabs(lastObj))+1.0e-4
4163      &&givenDuals==NULL&&firstFree_<0) {
4164    int maxFactor = factorization_->maximumPivots();
4165    if (maxFactor>10) {
4166      if (forceFactorization_<0)
4167        forceFactorization_= maxFactor;
4168      forceFactorization_ = CoinMax(1,(forceFactorization_>>1));
4169      //printf("Reducing factorization frequency\n");
4170    } 
4171  }
4172#endif
4173  // Allow matrices to be sorted etc
4174  int fake=-999; // signal sort
4175  matrix_->correctSequence(this,fake,fake);
4176  if (alphaAccuracy_>0.0)
4177      alphaAccuracy_=1.0;
4178}
4179/* While updateDualsInDual sees what effect is of flip
4180   this does actual flipping.
4181   If change >0.0 then value in array >0.0 => from lower to upper
4182*/
4183void 
4184ClpSimplexDual::flipBounds(CoinIndexedVector * rowArray,
4185                  CoinIndexedVector * columnArray,
4186                  double change)
4187{
4188  int number;
4189  int * which;
4190 
4191  int iSection;
4192
4193  for (iSection=0;iSection<2;iSection++) {
4194    int i;
4195    double * solution = solutionRegion(iSection);
4196    double * lower = lowerRegion(iSection);
4197    double * upper = upperRegion(iSection);
4198    int addSequence;
4199    if (!iSection) {
4200      number = rowArray->getNumElements();
4201      which = rowArray->getIndices();
4202      addSequence = numberColumns_;
4203    } else {
4204      number = columnArray->getNumElements();
4205      which = columnArray->getIndices();
4206      addSequence = 0;
4207    }
4208   
4209    for (i=0;i<number;i++) {
4210      int iSequence = which[i];
4211      Status status = getStatus(iSequence+addSequence);
4212
4213      switch(status) {
4214
4215      case basic:
4216      case isFree:
4217      case superBasic:
4218      case ClpSimplex::isFixed:
4219        break;
4220      case atUpperBound:
4221        // to lower bound
4222        setStatus(iSequence+addSequence,atLowerBound);
4223        solution[iSequence] = lower[iSequence];
4224        break;
4225      case atLowerBound:
4226        // to upper bound
4227        setStatus(iSequence+addSequence,atUpperBound);
4228        solution[iSequence] = upper[iSequence];
4229        break;
4230      }
4231    }
4232  }
4233  rowArray->setNumElements(0);
4234  columnArray->setNumElements(0);
4235}
4236// Restores bound to original bound
4237void 
4238ClpSimplexDual::originalBound( int iSequence)
4239{
4240  if (getFakeBound(iSequence)!=noFake) {
4241    numberFake_--;;
4242    setFakeBound(iSequence,noFake);
4243    if (auxiliaryModel_) {
4244      // just copy back
4245      lower_[iSequence]=auxiliaryModel_->lowerRegion()[iSequence+numberRows_+numberColumns_];
4246      upper_[iSequence]=auxiliaryModel_->upperRegion()[iSequence+numberRows_+numberColumns_];
4247      return;
4248    }
4249    if (iSequence>=numberColumns_) {
4250      // rows
4251      int iRow = iSequence-numberColumns_;
4252      rowLowerWork_[iRow]=rowLower_[iRow];
4253      rowUpperWork_[iRow]=rowUpper_[iRow];
4254      if (rowScale_) {
4255        if (rowLowerWork_[iRow]>-1.0e50)
4256          rowLowerWork_[iRow] *= rowScale_[iRow]*rhsScale_;
4257        if (rowUpperWork_[iRow]<1.0e50)
4258          rowUpperWork_[iRow] *= rowScale_[iRow]*rhsScale_;
4259      } else if (rhsScale_!=1.0) {
4260        if (rowLowerWork_[iRow]>-1.0e50)
4261          rowLowerWork_[iRow] *= rhsScale_;
4262        if (rowUpperWork_[iRow]<1.0e50)
4263          rowUpperWork_[iRow] *= rhsScale_;
4264      }
4265    } else {
4266      // columns
4267      columnLowerWork_[iSequence]=columnLower_[iSequence];
4268      columnUpperWork_[iSequence]=columnUpper_[iSequence];
4269      if (rowScale_) {
4270        double multiplier = 1.0/columnScale_[iSequence];
4271        if (columnLowerWork_[iSequence]>-1.0e50)
4272          columnLowerWork_[iSequence] *= multiplier*rhsScale_;
4273        if (columnUpperWork_[iSequence]<1.0e50)
4274          columnUpperWork_[iSequence] *= multiplier*rhsScale_;
4275      } else if (rhsScale_!=1.0) {
4276        if (columnLowerWork_[iSequence]>-1.0e50)
4277          columnLowerWork_[iSequence] *= rhsScale_;
4278        if (columnUpperWork_[iSequence]<1.0e50)
4279          columnUpperWork_[iSequence] *= rhsScale_;
4280      }
4281    }
4282  }
4283}
4284/* As changeBounds but just changes new bounds for a single variable.
4285   Returns true if change */
4286bool 
4287ClpSimplexDual::changeBound( int iSequence)
4288{
4289  // old values
4290  double oldLower=lower_[iSequence];
4291  double oldUpper=upper_[iSequence];
4292  double value=solution_[iSequence];
4293  bool modified=false;
4294  originalBound(iSequence);
4295  // original values
4296  double lowerValue=lower_[iSequence];
4297  double upperValue=upper_[iSequence];
4298  // back to altered values
4299  lower_[iSequence] = oldLower;
4300  upper_[iSequence] = oldUpper;
4301  assert (getFakeBound(iSequence)==noFake);
4302  //if (getFakeBound(iSequence)!=noFake)
4303  //numberFake_--;;
4304  if (value==oldLower) {
4305    if (upperValue > oldLower + dualBound_) {
4306      upper_[iSequence]=oldLower+dualBound_;
4307      setFakeBound(iSequence,upperFake);
4308      modified=true;
4309      numberFake_++;
4310    }
4311  } else if (value==oldUpper) {
4312    if (lowerValue < oldUpper - dualBound_) {
4313      lower_[iSequence]=oldUpper-dualBound_;
4314      setFakeBound(iSequence,lowerFake);
4315      modified=true;
4316      numberFake_++;
4317    }
4318  } else {
4319    assert(value==oldLower||value==oldUpper);
4320  }
4321  return modified;
4322}
4323// Perturbs problem
4324void 
4325ClpSimplexDual::perturb()
4326{
4327  if (perturbation_>100)
4328    return; //perturbed already
4329  if (perturbation_==100)
4330    perturbation_=50; // treat as normal
4331  int savePerturbation = perturbation_;
4332  bool modifyRowCosts=false;
4333  // dual perturbation
4334  double perturbation=1.0e-20;
4335  // maximum fraction of cost to perturb
4336  double maximumFraction = 1.0e-5;
4337  double constantPerturbation = 100.0*dualTolerance_;
4338  int maxLength=0;
4339  int minLength=numberRows_;
4340  double averageCost = 0.0;
4341  // look at element range
4342  double smallestNegative;
4343  double largestNegative;
4344  double smallestPositive;
4345  double largestPositive;
4346  matrix_->rangeOfElements(smallestNegative, largestNegative,
4347                           smallestPositive, largestPositive);
4348  smallestPositive = CoinMin(fabs(smallestNegative),smallestPositive);
4349  largestPositive = CoinMax(fabs(largestNegative),largestPositive);
4350  double elementRatio = largestPositive/smallestPositive;
4351  int numberNonZero=0;
4352  if (!numberIterations_&&perturbation_>=50) {
4353    // See if we need to perturb
4354    double * sort = new double[numberColumns_];
4355    // Use objective BEFORE scaling
4356    const double * obj = objective();
4357    int i;
4358    for (i=0;i<numberColumns_;i++) {
4359      double value = fabs(obj[i]);
4360      sort[i]=value;
4361      averageCost += value;
4362      if (value)
4363        numberNonZero++;
4364    }
4365    if (numberNonZero)
4366      averageCost /= (double) numberNonZero;
4367    else
4368      averageCost = 1.0;
4369    std::sort(sort,sort+numberColumns_);
4370    int number=1;
4371    double last = sort[0];
4372    for (i=1;i<numberColumns_;i++) {
4373      if (last!=sort[i])
4374        number++;
4375      last=sort[i];
4376    }
4377    delete [] sort;
4378#if 0
4379    printf("nnz %d percent %d",number,(number*100)/numberColumns_);
4380    if (number*4>numberColumns_)
4381      printf(" - Would not perturb\n");
4382    else
4383      printf(" - Would perturb\n");
4384    //exit(0);
4385#endif
4386    //printf("ratio number diff costs %g, element ratio %g\n",((double)number)/((double) numberColumns_),
4387    //                                                                elementRatio);
4388    //number=0;
4389    if (number*4>numberColumns_||elementRatio>1.0e12) {
4390      perturbation_=100;
4391      return; // good enough
4392    }
4393  }
4394  int iColumn;
4395  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4396    if (columnLowerWork_[iColumn]<columnUpperWork_[iColumn]) {
4397      int length = matrix_->getVectorLength(iColumn);
4398      if (length>2) {
4399        maxLength = CoinMax(maxLength,length);
4400        minLength = CoinMin(minLength,length);
4401      }
4402    }
4403  }
4404  // If > 70 then do rows
4405  if (perturbation_>=70) {
4406    modifyRowCosts=true;
4407    perturbation_ -= 20;
4408    printf("Row costs modified, ");
4409  }
4410  bool uniformChange=false;
4411  if (perturbation_>50) {
4412    // Experiment
4413    // maximumFraction could be 1.0e-10 to 1.0
4414    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};
4415    maximumFraction = m[CoinMin(perturbation_-51,10)];
4416  }
4417  int iRow;
4418  double smallestNonZero=1.0e100;
4419  numberNonZero=0;
4420  if (perturbation_>=50) {
4421    perturbation = 1.0e-8;
4422    bool allSame=true;
4423    double lastValue=0.0;
4424    for (iRow=0;iRow<numberRows_;iRow++) {
4425      double lo = rowLowerWork_[iRow];
4426      double up = rowUpperWork_[iRow];
4427      if (lo<up) {
4428        double value = fabs(rowObjectiveWork_[iRow]);
4429        perturbation = CoinMax(perturbation,value);
4430        if (value) {
4431          modifyRowCosts=true;
4432          smallestNonZero = CoinMin(smallestNonZero,value);
4433        }
4434      } 
4435      if (lo&&lo>-1.0e10) {
4436        numberNonZero++;
4437        lo=fabs(lo);
4438        if (!lastValue) 
4439          lastValue=lo;
4440        else if (fabs(lo-lastValue)>1.0e-7)
4441          allSame=false;
4442      }
4443      if (up&&up<1.0e10) {
4444        numberNonZero++;
4445        up=fabs(up);
4446        if (!lastValue) 
4447          lastValue=up;
4448        else if (fabs(up-lastValue)>1.0e-7)
4449          allSame=false;
4450      }
4451    }
4452    double lastValue2=0.0;
4453    for (iColumn=0;iColumn<numberColumns_;iColumn++) { 
4454      double lo = columnLowerWork_[iColumn];
4455      double up = columnUpperWork_[iColumn];
4456      if (lo<up) {
4457        double value = 
4458          fabs(objectiveWork_[iColumn]);
4459        perturbation = CoinMax(perturbation,value);
4460        if (value) {
4461          smallestNonZero = CoinMin(smallestNonZero,value);
4462        }
4463      }
4464      if (lo&&lo>-1.0e10) {
4465        //numberNonZero++;
4466        lo=fabs(lo);
4467        if (!lastValue2) 
4468          lastValue2=lo;
4469        else if (fabs(lo-lastValue2)>1.0e-7)
4470          allSame=false;
4471      }
4472      if (up&&up<1.0e10) {
4473        //numberNonZero++;
4474        up=fabs(up);
4475        if (!lastValue2) 
4476          lastValue2=up;
4477        else if (fabs(up-lastValue2)>1.0e-7)
4478          allSame=false;
4479      }
4480    }
4481    if (allSame) {
4482      // Check elements
4483      double smallestNegative;
4484      double largestNegative;
4485      double smallestPositive;
4486      double largestPositive;
4487      matrix_->rangeOfElements(smallestNegative,largestNegative,
4488                               smallestPositive,largestPositive);
4489      if (smallestNegative==largestNegative&&
4490          smallestPositive==largestPositive) {
4491        // Really hit perturbation
4492        double adjust = CoinMin(100.0*maximumFraction,1.0e-3*CoinMax(lastValue,lastValue2));
4493        maximumFraction = CoinMax(adjust,maximumFraction);
4494      }
4495    }
4496    perturbation = CoinMin(perturbation,smallestNonZero/maximumFraction);
4497  } else {
4498    // user is in charge
4499    maximumFraction = 1.0e-1;
4500    // but some experiments
4501    if (perturbation_<=-900) {
4502      modifyRowCosts=true;
4503      perturbation_ += 1000;
4504      printf("Row costs modified, ");
4505    }
4506    if (perturbation_<=-10) {
4507      perturbation_ += 10; 
4508      maximumFraction = 1.0;
4509      if ((-perturbation_)%100>=10) {
4510        uniformChange=true;
4511        perturbation_+=20;
4512      }
4513      while (perturbation_<-10) {
4514        perturbation_ += 100;
4515        maximumFraction *= 1.0e-1;
4516      }
4517    }
4518    perturbation = pow(10.0,perturbation_);
4519  }
4520  double largestZero=0.0;
4521  double largest=0.0;
4522  double largestPerCent=0.0;
4523  // modify costs
4524  bool printOut=(handler_->logLevel()==63);
4525  printOut=false;
4526  if (modifyRowCosts) {
4527    for (iRow=0;iRow<numberRows_;iRow++) {
4528      if (rowLowerWork_[iRow]<rowUpperWork_[iRow]) {
4529        double value = perturbation;
4530        double currentValue = rowObjectiveWork_[iRow];
4531        value = CoinMin(value,maximumFraction*(fabs(currentValue)+1.0e-1*perturbation+1.0e-3));
4532        if (rowLowerWork_[iRow]>-largeValue_) {
4533          if (fabs(rowLowerWork_[iRow])<fabs(rowUpperWork_[iRow])) 
4534            value *= CoinDrand48();
4535          else
4536            value *= -CoinDrand48();
4537        } else if (rowUpperWork_[iRow]<largeValue_) {
4538          value *= -CoinDrand48();
4539        } else {
4540          value=0.0;
4541        }
4542        if (currentValue) {
4543          largest = CoinMax(largest,fabs(value));
4544          if (fabs(value)>fabs(currentValue)*largestPerCent)
4545            largestPerCent=fabs(value/currentValue);
4546        } else {
4547          largestZero=CoinMax(largestZero,fabs(value));
4548        }
4549        if (printOut)
4550          printf("row %d cost %g change %g\n",iRow,rowObjectiveWork_[iRow],value);
4551        rowObjectiveWork_[iRow] += value;
4552      }
4553    }
4554  }
4555  // 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};
4556  // 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};
4557  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};
4558  //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};
4559  double extraWeight=10.0;
4560  // Scale back if wanted
4561  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};
4562  if (constantPerturbation<99.0*dualTolerance_) {
4563    perturbation *= 0.1;
4564    extraWeight=0.5;
4565    memcpy(weight,weight2,sizeof(weight2));
4566  }
4567  // adjust weights if all columns long
4568  double factor=1.0;
4569  if (maxLength) {
4570    factor = 3.0/(double) minLength;
4571  }
4572  // Make variables with more elements more expensive
4573  const double m1 = 0.5;
4574  double smallestAllowed = CoinMin(1.0e-2*dualTolerance_,maximumFraction);
4575  //double largestAllowed = CoinMax(1.0e3*dualTolerance_,maximumFraction*10.0*averageCost);
4576  double largestAllowed = CoinMax(1.0e3*dualTolerance_,maximumFraction*averageCost);
4577  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4578    if (columnLowerWork_[iColumn]<columnUpperWork_[iColumn]&&getStatus(iColumn)!=basic) {
4579      double value = perturbation;
4580      double currentValue = objectiveWork_[iColumn];
4581      value = CoinMin(value,constantPerturbation+maximumFraction*(fabs(currentValue)+1.0e-1*perturbation+1.0e-8));
4582      //value = CoinMin(value,constantPerturbation;+maximumFraction*fabs(currentValue));
4583      double value2 = constantPerturbation+1.0e-1*smallestNonZero;
4584      if (uniformChange) {
4585        value = maximumFraction;
4586        value2=maximumFraction;
4587      }
4588      if (columnLowerWork_[iColumn]>-largeValue_) {
4589        if (fabs(columnLowerWork_[iColumn])<
4590            fabs(columnUpperWork_[iColumn])) {
4591          value *= (1.0-m1 +m1*CoinDrand48());
4592          value2 *= (1.0-m1 +m1*CoinDrand48());
4593        } else {
4594          value *= -(1.0-m1+m1*CoinDrand48());
4595          value2 *= -(1.0-m1+m1*CoinDrand48());
4596        }
4597      } else if (columnUpperWork_[iColumn]<largeValue_) {
4598        value *= -(1.0-m1+m1*CoinDrand48());
4599        value2 *= -(1.0-m1+m1*CoinDrand48());
4600      } else {
4601        value=0.0;
4602      }
4603      if (value) {
4604        int length = matrix_->getVectorLength(iColumn);
4605        if (length>3) {
4606          length = (int) ((double) length * factor);
4607          length = CoinMax(3,length);
4608        }
4609        double multiplier;
4610#if 1
4611        if (length<10)
4612          multiplier=weight[length];
4613        else
4614          multiplier = weight[10];
4615#else
4616        if (length<10)
4617          multiplier=weight[length];
4618        else
4619          multiplier = weight[10]+extraWeight*(length-10);
4620        multiplier *= 0.5;
4621#endif
4622        value *= multiplier;
4623        value = CoinMin(value,value2);
4624        if (savePerturbation<50||savePerturbation>60) {
4625          if (fabs(value)<=dualTolerance_)
4626            value=0.0;
4627        } else if (value) {
4628          // get in range
4629          if (fabs(value)<=smallestAllowed) {
4630            value *= 10.0;
4631            while (fabs(value)<=smallestAllowed) 
4632              value *= 10.0;
4633          } else if (fabs(value)>largestAllowed) {
4634            value *= 0.1;
4635            while (fabs(value)>largestAllowed) 
4636              value *= 0.1;
4637          }
4638        }
4639        if (currentValue) {
4640          largest = CoinMax(largest,fabs(value));
4641          if (fabs(value)>fabs(currentValue)*largestPerCent)
4642            largestPerCent=fabs(value/currentValue);
4643        } else {
4644          largestZero=CoinMax(largestZero,fabs(value));
4645        }
4646        // but negative if at ub
4647        if (getStatus(iColumn)==atUpperBound)
4648          value = -value;
4649        if (printOut)
4650          printf("col %d cost %g change %g\n",iColumn,objectiveWork_[iColumn],value);
4651        objectiveWork_[iColumn] += value;
4652      }
4653    }
4654  }
4655  handler_->message(CLP_SIMPLEX_PERTURB,messages_)
4656    <<100.0*maximumFraction<<perturbation<<largest<<100.0*largestPerCent<<largestZero
4657    <<CoinMessageEol;
4658  // and zero changes
4659  //int nTotal = numberRows_+numberColumns_;
4660  //CoinZeroN(cost_+nTotal,nTotal);
4661  // say perturbed
4662  perturbation_=101;
4663
4664}
4665/* For strong branching.  On input lower and upper are new bounds
4666   while on output they are change in objective function values
4667   (>1.0e50 infeasible).
4668   Return code is 0 if nothing interesting, -1 if infeasible both
4669   ways and +1 if infeasible one way (check values to see which one(s))
4670   Returns -2 if bad factorization
4671*/
4672int ClpSimplexDual::strongBranching(int numberVariables,const int * variables,
4673                                    double * newLower, double * newUpper,
4674                                    double ** outputSolution,
4675                                    int * outputStatus, int * outputIterations,
4676                                    bool stopOnFirstInfeasible,
4677                                    bool alwaysFinish,
4678                                    int startFinishOptions)
4679{
4680  int i;
4681  int returnCode=0;
4682  double saveObjectiveValue = objectiveValue_;
4683  algorithm_ = -1;
4684
4685  //scaling(false);
4686
4687  // put in standard form (and make row copy)
4688  // create modifiable copies of model rim and do optional scaling
4689  createRim(7+8+16+32,true,startFinishOptions);
4690
4691  // change newLower and newUpper if scaled
4692
4693  // Do initial factorization
4694  // and set certain stuff
4695  // We can either set increasing rows so ...IsBasic gives pivot row
4696  // or we can just increment iBasic one by one
4697  // for now let ...iBasic give pivot row
4698  int useFactorization=false;
4699  if ((startFinishOptions&2)!=0&&(whatsChanged_&(2+512))==2+512)
4700    useFactorization=true; // Keep factorization if possible
4701  // switch off factorization if bad
4702  if (pivotVariable_[0]<0)
4703    useFactorization=false;
4704  if (!useFactorization||factorization_->numberRows()!=numberRows_) {
4705    useFactorization = false;
4706    factorization_->increasingRows(2);
4707    // row activities have negative sign
4708    factorization_->slackValue(-1.0);
4709    factorization_->zeroTolerance(1.0e-13);
4710
4711    int factorizationStatus = internalFactorize(0);
4712    if (factorizationStatus<0) {
4713      // some error
4714      // we should either debug or ignore
4715#ifndef NDEBUG
4716      printf("***** ClpDual strong branching factorization error - debug\n");
4717#endif
4718      return -2;
4719    } else if (factorizationStatus&&factorizationStatus<=numberRows_) {
4720      handler_->message(CLP_SINGULARITIES,messages_)
4721        <<factorizationStatus
4722        <<CoinMessageEol;
4723    }
4724  }
4725  // save stuff
4726  ClpFactorization saveFactorization(*factorization_);
4727  // save basis and solution
4728  double * saveSolution = new double[numberRows_+numberColumns_];
4729  CoinMemcpyN(solution_,
4730         numberRows_+numberColumns_,saveSolution);
4731  unsigned char * saveStatus = 
4732    new unsigned char [numberRows_+numberColumns_];
4733  CoinMemcpyN(status_,numberColumns_+numberRows_,saveStatus);
4734  // save bounds as createRim makes clean copies
4735  double * saveLower = new double[numberRows_+numberColumns_];
4736  CoinMemcpyN(lower_,
4737         numberRows_+numberColumns_,saveLower);
4738  double * saveUpper = new double[numberRows_+numberColumns_];
4739  CoinMemcpyN(upper_,
4740         numberRows_+numberColumns_,saveUpper);
4741  double * saveObjective = new double[numberRows_+numberColumns_];
4742  CoinMemcpyN(cost_,
4743         numberRows_+numberColumns_,saveObjective);
4744  int * savePivot = new int [numberRows_];
4745  CoinMemcpyN(pivotVariable_, numberRows_,savePivot);
4746  // need to save/restore weights.
4747
4748  int iSolution = 0;
4749  for (i=0;i<numberVariables;i++) {
4750    int iColumn = variables[i];
4751    double objectiveChange;
4752    double saveBound;
4753   
4754    // try down
4755   
4756    saveBound = columnUpper_[iColumn];
4757    // external view - in case really getting optimal
4758    columnUpper_[iColumn] = newUpper[i];
4759    if (scalingFlag_<=0) 
4760      upper_[iColumn] = newUpper[i]*rhsScale_;
4761    else 
4762      upper_[iColumn] = (newUpper[i]/columnScale_[iColumn])*rhsScale_; // scale
4763    // Start of fast iterations
4764    int status = fastDual(alwaysFinish);
4765    CoinAssert (problemStatus_||objectiveValue_<1.0e50);
4766    // make sure plausible
4767    double obj = CoinMax(objectiveValue_,saveObjectiveValue);
4768    if (status&&problemStatus_!=3) {
4769      // not finished - might be optimal
4770      checkPrimalSolution(rowActivityWork_,columnActivityWork_);
4771      double limit = 0.0;
4772      getDblParam(ClpDualObjectiveLimit, limit);
4773      if (!numberPrimalInfeasibilities_&&obj<limit) { 
4774        problemStatus_=0;
4775      } 
4776      status=problemStatus_;
4777    }
4778    if (status||(problemStatus_==0&&!isDualObjectiveLimitReached())) {
4779      objectiveChange = obj-saveObjectiveValue;
4780    } else {
4781      objectiveChange = 1.0e100;
4782      status=1;
4783    }
4784    if (problemStatus_==3)
4785      status=2;
4786
4787    if (scalingFlag_<=0) {
4788      CoinMemcpyN(solution_,numberColumns_,outputSolution[iSolution]);
4789    } else {
4790      int j;
4791      double * sol = outputSolution[iSolution];
4792      for (j=0;j<numberColumns_;j++) 
4793        sol[j] = solution_[j]*columnScale_[j];
4794    }
4795    outputStatus[iSolution]=status;
4796    outputIterations[iSolution]=numberIterations_;
4797    iSolution++;
4798    // restore
4799    CoinMemcpyN(saveSolution,
4800            numberRows_+numberColumns_,solution_);
4801    CoinMemcpyN(saveStatus,numberColumns_+numberRows_,status_);
4802    CoinMemcpyN(saveLower,
4803            numberRows_+numberColumns_,lower_);
4804    CoinMemcpyN(saveUpper,
4805            numberRows_+numberColumns_,upper_);
4806    CoinMemcpyN(saveObjective,
4807            numberRows_+numberColumns_,cost_);
4808    columnUpper_[iColumn] = saveBound;
4809    CoinMemcpyN(savePivot, numberRows_,pivotVariable_);
4810    delete factorization_;
4811    factorization_ = new ClpFactorization(saveFactorization);
4812
4813    newUpper[i]=objectiveChange;
4814#ifdef CLP_DEBUG
4815    printf("down on %d costs %g\n",iColumn,objectiveChange);
4816#endif
4817         
4818    // try up
4819   
4820    saveBound = columnLower_[iColumn];
4821    // external view - in case really getting optimal
4822    columnLower_[iColumn] = newLower[i];
4823    if (scalingFlag_<=0) 
4824      lower_[iColumn] = newLower[i]*rhsScale_;
4825    else 
4826      lower_[iColumn] = (newLower[i]/columnScale_[iColumn])*rhsScale_; // scale
4827    // Start of fast iterations
4828    status = fastDual(alwaysFinish);
4829    // make sure plausible
4830    obj = CoinMax(objectiveValue_,saveObjectiveValue);
4831    if (status&&problemStatus_!=3) {
4832      // not finished - might be optimal
4833      checkPrimalSolution(rowActivityWork_,columnActivityWork_);
4834      double limit = 0.0;
4835      getDblParam(ClpDualObjectiveLimit, limit);
4836      if (!numberPrimalInfeasibilities_&&obj< limit) { 
4837        problemStatus_=0;
4838      } 
4839      status=problemStatus_;
4840    }
4841    if (status||(problemStatus_==0&&!isDualObjectiveLimitReached())) {
4842      objectiveChange = obj-saveObjectiveValue;
4843    } else {
4844      objectiveChange = 1.0e100;
4845      status=1;
4846    }
4847    if (problemStatus_==3)
4848      status=2;
4849    if (scalingFlag_<=0) {
4850      CoinMemcpyN(solution_,numberColumns_,outputSolution[iSolution]);
4851    } else {
4852      int j;
4853      double * sol = outputSolution[iSolution];
4854      for (j=0;j<numberColumns_;j++) 
4855        sol[j] = solution_[j]*columnScale_[j];
4856    }
4857    outputStatus[iSolution]=status;
4858    outputIterations[iSolution]=numberIterations_;
4859    iSolution++;
4860
4861    // restore
4862    CoinMemcpyN(saveSolution,
4863            numberRows_+numberColumns_,solution_);
4864    CoinMemcpyN(saveStatus,numberColumns_+numberRows_,status_);
4865    CoinMemcpyN(saveLower,
4866            numberRows_+numberColumns_,lower_);
4867    CoinMemcpyN(saveUpper,
4868            numberRows_+numberColumns_,upper_);
4869    CoinMemcpyN(saveObjective,
4870            numberRows_+numberColumns_,cost_);
4871    columnLower_[iColumn] = saveBound;
4872    CoinMemcpyN(savePivot, numberRows_,pivotVariable_);
4873    delete factorization_;
4874    factorization_ = new ClpFactorization(saveFactorization);
4875
4876    newLower[i]=objectiveChange;
4877#ifdef CLP_DEBUG
4878    printf("up on %d costs %g\n",iColumn,objectiveChange);
4879#endif
4880         
4881    /* Possibilities are:
4882       Both sides feasible - store
4883       Neither side feasible - set objective high and exit if desired
4884       One side feasible - change bounds and resolve
4885    */
4886    if (newUpper[i]<1.0e100) {
4887      if(newLower[i]<1.0e100) {
4888        // feasible - no action
4889      } else {
4890        // up feasible, down infeasible
4891        returnCode=1;
4892        if (stopOnFirstInfeasible)
4893          break;
4894      }
4895    } else {
4896      if(newLower[i]<1.0e100) {
4897        // down feasible, up infeasible
4898        returnCode=1;
4899        if (stopOnFirstInfeasible)
4900          break;
4901      } else {
4902        // neither side feasible
4903        returnCode=-1;
4904        break;
4905      }
4906    }
4907  }
4908  delete [] saveSolution;
4909  delete [] saveLower;
4910  delete [] saveUpper;
4911  delete [] saveObjective;
4912  delete [] saveStatus;
4913  delete [] savePivot;
4914  if ((startFinishOptions&1)==0) {
4915    deleteRim(1);
4916    whatsChanged_=0;
4917  } else {
4918    // Original factorization will have been put back by last loop
4919    //delete factorization_;
4920    //factorization_ = new ClpFactorization(saveFactorization);
4921    deleteRim(0);
4922    // mark all as current
4923    whatsChanged_ = 0xffff;
4924  }
4925  objectiveValue_ = saveObjectiveValue;
4926  return returnCode;
4927}
4928// treat no pivot as finished (unless interesting)
4929int ClpSimplexDual::fastDual(bool alwaysFinish)
4930{
4931  algorithm_ = -1;
4932  secondaryStatus_=0;
4933  // Say in fast dual
4934  specialOptions_ |= 16384;
4935  //handler_->setLogLevel(63);
4936  // save data
4937  ClpDataSave data = saveData();
4938  dualTolerance_=dblParam_[ClpDualTolerance];
4939  primalTolerance_=dblParam_[ClpPrimalTolerance];
4940
4941  // save dual bound
4942  double saveDualBound = dualBound_;
4943
4944  if (alphaAccuracy_!=-1.0)
4945    alphaAccuracy_ = 1.0;
4946  double objectiveChange;
4947  // for dual we will change bounds using dualBound_
4948  // for this we need clean basis so it is after factorize
4949  gutsOfSolution(NULL,NULL);
4950  numberFake_ =0; // Number of variables at fake bounds
4951  numberChanged_ =0; // Number of variables with changed costs
4952  changeBounds(true,NULL,objectiveChange);
4953
4954  problemStatus_ = -1;
4955  numberIterations_=0;
4956  factorization_->sparseThreshold(0);
4957  factorization_->goSparse();
4958
4959  int lastCleaned=0; // last time objective or bounds cleaned up
4960
4961  // number of times we have declared optimality
4962  numberTimesOptimal_=0;
4963
4964  // This says whether to restore things etc
4965  int factorType=0;
4966  /*
4967    Status of problem:
4968    0 - optimal
4969    1 - infeasible
4970    2 - unbounded
4971    -1 - iterating
4972    -2 - factorization wanted
4973    -3 - redo checking without factorization
4974    -4 - looks infeasible
4975
4976    BUT also from whileIterating return code is:
4977
4978   -1 iterations etc
4979   -2 inaccuracy
4980   -3 slight inaccuracy (and done iterations)
4981   +0 looks optimal (might be unbounded - but we will investigate)
4982   +1 looks infeasible
4983   +3 max iterations
4984
4985  */
4986
4987  int returnCode = 0;
4988
4989  int iRow,iColumn;
4990  while (problemStatus_<0) {
4991    // clear
4992    for (iRow=0;iRow<4;iRow++) {
4993      rowArray_[iRow]->clear();
4994    }   
4995   
4996    for (iColumn=0;iColumn<2;iColumn++) {
4997      columnArray_[iColumn]->clear();
4998    }   
4999
5000    // give matrix (and model costs and bounds a chance to be
5001    // refreshed (normally null)
5002    matrix_->refresh(this);
5003    // may factorize, checks if problem finished
5004    // should be able to speed this up on first time
5005    statusOfProblemInDual(lastCleaned,factorType,NULL,data,0);
5006
5007    // Say good factorization
5008    factorType=1;
5009
5010    // Do iterations
5011    if (problemStatus_<0) {
5012      double * givenPi=NULL;
5013      returnCode = whileIterating(givenPi,0);
5014      if ((!alwaysFinish&&returnCode<0)||returnCode==3) {
5015        if (returnCode!=3)
5016          assert (problemStatus_<0);
5017        returnCode=1;
5018        problemStatus_ = 3; 
5019        // can't say anything interesting - might as well return
5020#ifdef CLP_DEBUG
5021        printf("returning from fastDual after %d iterations with code %d\n",
5022               numberIterations_,returnCode);
5023#endif
5024        break;
5025      }
5026      returnCode=0;
5027    }
5028  }
5029
5030  // clear
5031  for (iRow=0;iRow<4;iRow++) {
5032    rowArray_[iRow]->clear();
5033  }   
5034 
5035  for (iColumn=0;iColumn<2;iColumn++) {
5036    columnArray_[iColumn]->clear();
5037  }   
5038  // Say not in fast dual
5039  specialOptions_ &= ~16384;
5040  assert(!numberFake_||((specialOptions_&(2048|4096))!=0&&dualBound_>1.0e8)
5041         ||returnCode||problemStatus_); // all bounds should be okay
5042  // Restore any saved stuff
5043  restoreData(data);
5044  dualBound_ = saveDualBound;
5045  return returnCode;
5046}
5047// This does first part of StrongBranching
5048ClpFactorization * 
5049ClpSimplexDual::setupForStrongBranching(char * arrays, int numberRows, int numberColumns)
5050{
5051  algorithm_ = -1;
5052  // put in standard form (and make row copy)
5053  // create modifiable copies of model rim and do optional scaling
5054  int startFinishOptions;
5055  /*  COIN_CLP_VETTED
5056      Looks safe for Cbc
5057  */
5058  bool safeOption = ((specialOptions_&COIN_CBC_USING_CLP)!=0);
5059  if((specialOptions_&4096)==0||(auxiliaryModel_&&!safeOption)) {
5060    startFinishOptions=0;
5061  } else {
5062    startFinishOptions=1+2+4;
5063  }
5064  createRim(7+8+16+32,true,startFinishOptions);
5065  // Do initial factorization
5066  // and set certain stuff
5067  // We can either set increasing rows so ...IsBasic gives pivot row
5068  // or we can just increment iBasic one by one
5069  // for now let ...iBasic give pivot row
5070  bool useFactorization=false;
5071  if ((startFinishOptions&2)!=0&&(whatsChanged_&(2+512))==2+512) {
5072    useFactorization=true; // Keep factorization if possible
5073    // switch off factorization if bad
5074    if (pivotVariable_[0]<0||factorization_->numberRows()!=numberRows_) 
5075      useFactorization=false;
5076  }
5077  if (!useFactorization) {
5078    factorization_->increasingRows(2);
5079    // row activities have negative sign
5080    factorization_->slackValue(-1.0);
5081    factorization_->zeroTolerance(1.0e-13);
5082
5083    int factorizationStatus = internalFactorize(0);
5084    if (factorizationStatus<0) {
5085      // some error
5086      // we should either debug or ignore
5087#ifndef NDEBUG
5088      printf("***** ClpDual strong branching factorization error - debug\n");
5089#endif
5090    } else if (factorizationStatus&&factorizationStatus<=numberRows_) {
5091      handler_->message(CLP_SINGULARITIES,messages_)
5092        <<factorizationStatus
5093        <<CoinMessageEol;
5094    }
5095  }
5096  double * arrayD = (double *) arrays;
5097  arrayD[0]=objectiveValue()*optimizationDirection_;
5098  double * saveSolution = arrayD+1;
5099  double * saveLower = saveSolution + (numberRows+numberColumns);
5100  double * saveUpper = saveLower + (numberRows+numberColumns);
5101  double * saveObjective = saveUpper + (numberRows+numberColumns);
5102  double * saveLowerOriginal = saveObjective + (numberRows+numberColumns);
5103  double * saveUpperOriginal = saveLowerOriginal + numberColumns;
5104  arrayD = saveUpperOriginal + numberColumns;
5105  int * savePivot = (int *) arrayD;
5106  int * whichRow = savePivot+numberRows;
5107  int * whichColumn = whichRow + 3*numberRows;
5108  int * arrayI = whichColumn + 2*numberColumns;
5109  unsigned char * saveStatus = (unsigned char *) (arrayI+1);
5110  // save stuff
5111  // save basis and solution
5112  CoinMemcpyN(solution_,
5113          numberRows_+numberColumns_,saveSolution);
5114  CoinMemcpyN(status_,numberColumns_+numberRows_,saveStatus);
5115  CoinMemcpyN(lower_,
5116          numberRows_+numberColumns_,saveLower);
5117  CoinMemcpyN(upper_,
5118          numberRows_+numberColumns_,saveUpper);
5119  CoinMemcpyN(cost_,
5120          numberRows_+numberColumns_,saveObjective);
5121  CoinMemcpyN(pivotVariable_, numberRows_,savePivot);
5122  return new ClpFactorization(*factorization_);
5123}
5124// This cleans up after strong branching
5125void 
5126ClpSimplexDual::cleanupAfterStrongBranching()
5127{
5128  int startFinishOptions;
5129  /*  COIN_CLP_VETTED
5130      Looks safe for Cbc
5131  */
5132  bool safeOption = ((specialOptions_&COIN_CBC_USING_CLP)!=0);
5133  if((specialOptions_&4096)==0||(auxiliaryModel_&&!safeOption)) {
5134    startFinishOptions=0;
5135  } else {
5136    startFinishOptions=1+2+4;
5137  }
5138  if ((startFinishOptions&1)==0) {
5139    deleteRim(1);
5140    whatsChanged_=0;
5141  } else {
5142    // Original factorization will have been put back by last loop
5143    //delete factorization_;
5144    //factorization_ = new ClpFactorization(saveFactorization);
5145    deleteRim(0);
5146    // mark all as current
5147    whatsChanged_ = 0xffff;
5148  }
5149}
5150/* Checks number of variables at fake bounds.  This is used by fastDual
5151   so can exit gracefully before end */
5152int 
5153ClpSimplexDual::numberAtFakeBound()
5154{
5155  int iSequence;
5156  int numberFake=0;
5157 
5158  for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
5159    FakeBound bound = getFakeBound(iSequence);
5160    switch(getStatus(iSequence)) {
5161
5162    case basic:
5163      break;
5164    case isFree:
5165    case superBasic:
5166    case ClpSimplex::isFixed:
5167      //setFakeBound (iSequence, noFake);
5168      break;
5169    case atUpperBound:
5170      if (bound==upperFake||bound==bothFake)
5171        numberFake++;
5172      break;
5173    case atLowerBound:
5174      if (bound==lowerFake||bound==bothFake)
5175        numberFake++;
5176      break;
5177    }
5178  }
5179  numberFake_ = numberFake;
5180  return numberFake;
5181}
5182/* Pivot out a variable and choose an incoing one.  Assumes dual
5183   feasible - will not go through a reduced cost. 
5184   Returns step length in theta
5185   Returns ray in ray_ (or NULL if no pivot)
5186   Return codes as before but -1 means no acceptable pivot
5187*/
5188int 
5189ClpSimplexDual::pivotResult()
5190{
5191  abort();
5192  return 0;
5193}
5194/*
5195   Row array has row part of pivot row
5196   Column array has column part.
5197   This is used in dual values pass
5198*/
5199void
5200ClpSimplexDual::checkPossibleValuesMove(CoinIndexedVector * rowArray,
5201                                        CoinIndexedVector * columnArray,
5202                                        double acceptablePivot)
5203{
5204  double * work;
5205  int number;
5206  int * which;
5207  int iSection;
5208
5209  double tolerance = dualTolerance_*1.001;
5210
5211  double thetaDown = 1.0e31;
5212  double changeDown ;
5213  double thetaUp = 1.0e31;
5214  double bestAlphaDown = acceptablePivot*0.99999;
5215  double bestAlphaUp = acceptablePivot*0.99999;
5216  int sequenceDown =-1;
5217  int sequenceUp = sequenceOut_;
5218
5219  double djBasic = dj_[sequenceOut_];
5220  if (djBasic>0.0) {
5221    // basic at lower bound so directionOut_ 1 and -1 in pivot row
5222    // dj will go to zero on other way
5223    thetaUp = djBasic;
5224    changeDown = -lower_[sequenceOut_];
5225  } else {
5226    // basic at upper bound so directionOut_ -1 and 1 in pivot row
5227    // dj will go to zero on other way
5228    thetaUp = -djBasic;
5229    changeDown = upper_[sequenceOut_];
5230  }
5231  bestAlphaUp = 1.0;
5232  int addSequence;
5233
5234  double alphaUp=0.0;
5235  double alphaDown=0.0;
5236
5237  for (iSection=0;iSection<2;iSection++) {
5238
5239    int i;
5240    if (!iSection) {
5241      work = rowArray->denseVector();
5242      number = rowArray->getNumElements();
5243      which = rowArray->getIndices();
5244      addSequence = numberColumns_;
5245    } else {
5246      work = columnArray->denseVector();
5247      number = columnArray->getNumElements();
5248      which = columnArray->getIndices();
5249      addSequence = 0;
5250    }
5251   
5252    for (i=0;i<number;i++) {
5253      int iSequence = which[i];
5254      int iSequence2 = iSequence + addSequence;
5255      double alpha;
5256      double oldValue;
5257      double value;
5258
5259      switch(getStatus(iSequence2)) {
5260         
5261      case basic: 
5262        break;
5263      case ClpSimplex::isFixed:
5264        alpha = work[i];
5265        changeDown += alpha*upper_[iSequence2];
5266        break;
5267      case isFree:
5268      case superBasic:
5269        alpha = work[i];
5270        // dj must be effectively zero as dual feasible
5271        if (fabs(alpha)>bestAlphaUp) {
5272          thetaDown = 0.0;
5273          thetaUp = 0.0;
5274          bestAlphaDown = fabs(alpha);
5275          bestAlphaUp = bestAlphaDown;
5276          sequenceDown =iSequence2;
5277          sequenceUp = sequenceDown;
5278          alphaUp = alpha;
5279          alphaDown = alpha;
5280        }
5281        break;
5282      case atUpperBound:
5283        alpha = work[i];
5284        oldValue = dj_[iSequence2];
5285        changeDown += alpha*upper_[iSequence2];
5286        if (alpha>=acceptablePivot) {
5287          // might do other way
5288          value = oldValue+thetaUp*alpha;
5289          if (value>-tolerance) {
5290            if (value>tolerance||fabs(alpha)>bestAlphaUp) {
5291              thetaUp = -oldValue/alpha;
5292              bestAlphaUp = fabs(alpha);
5293              sequenceUp = iSequence2;
5294              alphaUp = alpha;
5295            }
5296          }
5297        } else if (alpha<=-acceptablePivot) {
5298          // might do this way
5299          value = oldValue-thetaDown*alpha;
5300          if (value>-tolerance) {
5301            if (value>tolerance||fabs(alpha)>bestAlphaDown) {
5302              thetaDown = oldValue/alpha;
5303              bestAlphaDown = fabs(alpha);
5304              sequenceDown = iSequence2;
5305              alphaDown = alpha;
5306            }
5307          }
5308        }
5309        break;
5310      case atLowerBound:
5311        alpha = work[i];
5312        oldValue = dj_[iSequence2];
5313        changeDown += alpha*lower_[iSequence2];
5314        if (alpha<=-acceptablePivot) {
5315          // might do other way
5316          value = oldValue+thetaUp*alpha;
5317          if (value<tolerance) {
5318            if (value<-tolerance||fabs(alpha)>bestAlphaUp) {
5319              thetaUp = -oldValue/alpha;
5320              bestAlphaUp = fabs(alpha);
5321              sequenceUp = iSequence2;
5322              alphaUp = alpha;
5323            }
5324          }
5325        } else if (alpha>=acceptablePivot) {
5326          // might do this way
5327          value = oldValue-thetaDown*alpha;
5328          if (value<tolerance) {
5329            if (value<-tolerance||fabs(alpha)>bestAlphaDown) {
5330              thetaDown = oldValue/alpha;
5331              bestAlphaDown = fabs(alpha);
5332              sequenceDown = iSequence2;
5333              alphaDown = alpha;
5334            }
5335          }
5336        }
5337        break;
5338      }
5339    }
5340  }
5341  thetaUp *= -1.0;
5342  double changeUp = -thetaUp * changeDown;
5343  changeDown = -thetaDown*changeDown;
5344  if (CoinMax(fabs(thetaDown),fabs(thetaUp))<1.0e-8) {
5345    // largest
5346    if (fabs(alphaDown)<fabs(alphaUp)) {
5347      sequenceDown =-1;
5348    }
5349  }
5350  // choose
5351  sequenceIn_=-1;
5352  if (changeDown>changeUp&&sequenceDown>=0) {
5353    theta_ = thetaDown;
5354    if (fabs(changeDown)<1.0e30)
5355      sequenceIn_ = sequenceDown;
5356    alpha_ = alphaDown;
5357#ifdef CLP_DEBUG
5358    if ((handler_->logLevel()&32))
5359      printf("predicted way - dirout %d, change %g,%g theta %g\n",
5360             directionOut_,changeDown,changeUp,theta_);
5361#endif
5362  } else {
5363    theta_ = thetaUp;
5364    if (fabs(changeUp)<1.0e30)
5365      sequenceIn_ = sequenceUp;
5366    alpha_ = alphaUp;
5367    if (sequenceIn_!=sequenceOut_) {
5368#ifdef CLP_DEBUG
5369      if ((handler_->logLevel()&32))
5370        printf("opposite way - dirout %d, change %g,%g theta %g\n",
5371               directionOut_,changeDown,changeUp,theta_);
5372#endif
5373    } else {
5374#ifdef CLP_DEBUG
5375      if ((handler_->logLevel()&32))
5376        printf("opposite way to zero dj - dirout %d, change %g,%g theta %g\n",
5377               directionOut_,changeDown,changeUp,theta_);
5378#endif
5379    }
5380  }
5381  if (sequenceIn_>=0) {
5382    lowerIn_ = lower_[sequenceIn_];
5383    upperIn_ = upper_[sequenceIn_];
5384    valueIn_ = solution_[sequenceIn_];
5385    dualIn_ = dj_[sequenceIn_];
5386
5387    if (alpha_<0.0) {
5388      // as if from upper bound
5389      directionIn_=-1;
5390      upperIn_=valueIn_;
5391    } else {
5392      // as if from lower bound
5393      directionIn_=1;
5394      lowerIn_=valueIn_;
5395    }
5396  }
5397}
5398/*
5399   Row array has row part of pivot row
5400   Column array has column part.
5401   This is used in cleanup
5402*/
5403void
5404ClpSimplexDual::checkPossibleCleanup(CoinIndexedVector * rowArray,
5405                                        CoinIndexedVector * columnArray,
5406                                        double acceptablePivot)
5407{
5408  double * work;
5409  int number;
5410  int * which;
5411  int iSection;
5412
5413  double tolerance = dualTolerance_*1.001;
5414
5415  double thetaDown = 1.0e31;
5416  double thetaUp = 1.0e31;
5417  double bestAlphaDown = acceptablePivot*10.0;
5418  double bestAlphaUp = acceptablePivot*10.0;
5419  int sequenceDown =-1;
5420  int sequenceUp = -1;
5421
5422  double djSlack = dj_[pivotRow_];
5423  if (getRowStatus(pivotRow_)==basic)
5424    djSlack=COIN_DBL_MAX;
5425  if (fabs(djSlack)<tolerance)
5426    djSlack=0.0;
5427  int addSequence;
5428
5429  double alphaUp=0.0;
5430  double alphaDown=0.0;
5431  for (iSection=0;iSection<2;iSection++) {
5432
5433    int i;
5434    if (!iSection) {
5435      work = rowArray->denseVector();
5436      number = rowArray->getNumElements();
5437      which = rowArray->getIndices();
5438      addSequence = numberColumns_;
5439    } else {
5440      work = columnArray->denseVector();
5441      number = columnArray->getNumElements();
5442      which = columnArray->getIndices();
5443      addSequence = 0;
5444    }
5445   
5446    for (i=0;i<number;i++) {
5447      int iSequence = which[i];
5448      int iSequence2 = iSequence + addSequence;
5449      double alpha;
5450      double oldValue;
5451      double value;
5452
5453      switch(getStatus(iSequence2)) {
5454         
5455      case basic: 
5456        break;
5457      case ClpSimplex::isFixed:
5458        alpha = work[i];
5459        if (addSequence) {
5460          printf("possible - pivot row %d this %d\n",pivotRow_,iSequence);
5461          oldValue = dj_[iSequence2];
5462          if (alpha<=-acceptablePivot) {
5463            // might do other way
5464            value = oldValue+thetaUp*alpha;
5465            if (value<tolerance) {
5466              if (value<-tolerance||fabs(alpha)>bestAlphaUp) {
5467                thetaUp = -oldValue/alpha;
5468                bestAlphaUp = fabs(alpha);
5469                sequenceUp = iSequence2;
5470                alphaUp = alpha;
5471              }
5472            }
5473          } else if (alpha>=acceptablePivot) {
5474            // might do this way
5475            value = oldValue-thetaDown*alpha;
5476            if (value<tolerance) {
5477              if (value<-tolerance||fabs(alpha)>bestAlphaDown) {
5478                thetaDown = oldValue/alpha;
5479                bestAlphaDown = fabs(alpha);
5480                sequenceDown = iSequence2;
5481                alphaDown = alpha;
5482              }
5483            }
5484          }
5485        }
5486        break;
5487      case isFree:
5488      case superBasic:
5489        alpha = work[i];
5490        // dj must be effectively zero as dual feasible
5491        if (fabs(alpha)>bestAlphaUp) {
5492          thetaDown = 0.0;
5493          thetaUp = 0.0;
5494          bestAlphaDown = fabs(alpha);
5495          bestAlphaUp = bestAlphaDown;
5496          sequenceDown =iSequence2;
5497          sequenceUp = sequenceDown;
5498          alphaUp = alpha;
5499          alphaDown = alpha;
5500        }
5501        break;
5502      case atUpperBound:
5503        alpha = work[i];
5504        oldValue = dj_[iSequence2];
5505        if (alpha>=acceptablePivot) {
5506          // might do other way
5507          value = oldValue+thetaUp*alpha;
5508          if (value>-tolerance) {
5509            if (value>tolerance||fabs(alpha)>bestAlphaUp) {
5510              thetaUp = -oldValue/alpha;
5511              bestAlphaUp = fabs(alpha);
5512              sequenceUp = iSequence2;
5513              alphaUp = alpha;
5514            }
5515          }
5516        } else if (alpha<=-acceptablePivot) {
5517          // might do this way
5518          value = oldValue-thetaDown*alpha;
5519          if (value>-tolerance) {
5520            if (value>tolerance||fabs(alpha)>bestAlphaDown) {
5521              thetaDown = oldValue/alpha;
5522              bestAlphaDown = fabs(alpha);
5523              sequenceDown = iSequence2;
5524              alphaDown = alpha;
5525            }
5526          }
5527        }
5528        break;
5529      case atLowerBound:
5530        alpha = work[i];
5531        oldValue = dj_[iSequence2];
5532        if (alpha<=-acceptablePivot) {
5533          // might do other way
5534          value = oldValue+thetaUp*alpha;
5535          if (value<tolerance) {
5536            if (value<-tolerance||fabs(alpha)>bestAlphaUp) {
5537              thetaUp = -oldValue/alpha;
5538              bestAlphaUp = fabs(alpha);
5539              sequenceUp = iSequence2;
5540              alphaUp = alpha;
5541            }
5542          }
5543        } else if (alpha>=acceptablePivot) {
5544          // might do this way
5545          value = oldValue-thetaDown*alpha;
5546          if (value<tolerance) {
5547            if (value<-tolerance||fabs(alpha)>bestAlphaDown) {
5548              thetaDown = oldValue/alpha;
5549              bestAlphaDown = fabs(alpha);
5550              sequenceDown = iSequence2;
5551              alphaDown = alpha;
5552            }
5553          }
5554        }
5555        break;
5556      }
5557    }
5558  }
5559  thetaUp *= -1.0;
5560  // largest
5561  if (bestAlphaDown<bestAlphaUp) 
5562    sequenceDown =-1;
5563  else
5564    sequenceUp=-1;
5565
5566  sequenceIn_=-1;
5567 
5568  if (sequenceDown>=0) {
5569    theta_ = thetaDown;
5570    sequenceIn_ = sequenceDown;
5571    alpha_ = alphaDown;
5572#ifdef CLP_DEBUG
5573    if ((handler_->logLevel()&32))
5574      printf("predicted way - dirout %d, theta %g\n",
5575             directionOut_,theta_);
5576#endif
5577  } else if (sequenceUp>=0) {
5578    theta_ = thetaUp;
5579    sequenceIn_ = sequenceUp;
5580    alpha_ = alphaUp;
5581#ifdef CLP_DEBUG
5582    if ((handler_->logLevel()&32))
5583      printf("opposite way - dirout %d,theta %g\n",
5584             directionOut_,theta_);
5585#endif
5586  }
5587  if (sequenceIn_>=0) {
5588    lowerIn_ = lower_[sequenceIn_];
5589    upperIn_ = upper_[sequenceIn_];
5590    valueIn_ = solution_[sequenceIn_];
5591    dualIn_ = dj_[sequenceIn_];
5592
5593    if (alpha_<0.0) {
5594      // as if from upper bound
5595      directionIn_=-1;
5596      upperIn_=valueIn_;
5597    } else {
5598      // as if from lower bound
5599      directionIn_=1;
5600      lowerIn_=valueIn_;
5601    }
5602  }
5603}
5604/*
5605   This sees if we can move duals in dual values pass.
5606   This is done before any pivoting
5607*/
5608void ClpSimplexDual::doEasyOnesInValuesPass(double * dj)
5609{
5610  // Get column copy
5611  CoinPackedMatrix * columnCopy = matrix();
5612  // Get a row copy in standard format
5613  CoinPackedMatrix copy;
5614  copy.reverseOrderedCopyOf(*columnCopy);
5615  // get matrix data pointers
5616  const int * column = copy.getIndices();
5617  const CoinBigIndex * rowStart = copy.getVectorStarts();
5618  const int * rowLength = copy.getVectorLengths(); 
5619  const double * elementByRow = copy.getElements();
5620  double tolerance = dualTolerance_*1.001;
5621
5622  int iRow;
5623#ifdef CLP_DEBUG
5624  {
5625    double value5=0.0;
5626    int i;
5627    for (i=0;i<numberRows_+numberColumns_;i++) {
5628      if (dj[i]<-1.0e-6)
5629        value5 += dj[i]*upper_[i];
5630      else if (dj[i] >1.0e-6)
5631        value5 += dj[i]*lower_[i];
5632    }
5633    printf("Values objective Value before %g\n",value5);
5634  }
5635#endif
5636  // for scaled row
5637  double * scaled=NULL;
5638  if (rowScale_)
5639    scaled = new double[numberColumns_];
5640  for (iRow=0;iRow<numberRows_;iRow++) {
5641
5642    int iSequence = iRow + numberColumns_;
5643    double djBasic = dj[iSequence];
5644    if (getRowStatus(iRow)==basic&&fabs(djBasic)>tolerance) {
5645
5646      double changeUp ;
5647      // always -1 in pivot row
5648      if (djBasic>0.0) {
5649        // basic at lower bound
5650        changeUp = -lower_[iSequence];
5651      } else {
5652        // basic at upper bound
5653        changeUp = upper_[iSequence];
5654      }
5655      bool canMove=true;
5656      int i;
5657      const double * thisElements = elementByRow + rowStart[iRow]; 
5658      const int * thisIndices = column+rowStart[iRow];
5659      if (rowScale_) {
5660        assert (!auxiliaryModel_);
5661        // scale row
5662        double scale = rowScale_[iRow];
5663        for (i = 0;i<rowLength[iRow];i++) {
5664          int iColumn = thisIndices[i];
5665          double alpha = thisElements[i];
5666          scaled[i] = scale*alpha*columnScale_[iColumn];
5667        }
5668        thisElements = scaled;
5669      }
5670      for (i = 0;i<rowLength[iRow];i++) {
5671        int iColumn = thisIndices[i];
5672        double alpha = thisElements[i];
5673        double oldValue = dj[iColumn];;
5674        double value;
5675
5676        switch(getStatus(iColumn)) {
5677         
5678        case basic:
5679          if (dj[iColumn]<-tolerance&&
5680              fabs(solution_[iColumn]-upper_[iColumn])<1.0e-8) {
5681            // at ub
5682            changeUp += alpha*upper_[iColumn];
5683            // might do other way
5684            value = oldValue+djBasic*alpha;
5685            if (value>tolerance) 
5686              canMove=false;
5687          } else if (dj[iColumn]>tolerance&&
5688              fabs(solution_[iColumn]-lower_[iColumn])<1.0e-8) {
5689            changeUp += alpha*lower_[iColumn];
5690            // might do other way
5691            value = oldValue+djBasic*alpha;
5692            if (value<-tolerance)
5693              canMove=false;
5694          } else {
5695            canMove=false;
5696          }
5697          break;
5698        case ClpSimplex::isFixed:
5699          changeUp += alpha*upper_[iColumn];
5700          break;
5701        case isFree:
5702        case superBasic:
5703          canMove=false;
5704        break;
5705      case atUpperBound:
5706        changeUp += alpha*upper_[iColumn];
5707        // might do other way
5708        value = oldValue+djBasic*alpha;
5709        if (value>tolerance) 
5710          canMove=false;
5711        break;
5712      case atLowerBound:
5713        changeUp += alpha*lower_[iColumn];
5714        // might do other way
5715        value = oldValue+djBasic*alpha;
5716        if (value<-tolerance)
5717          canMove=false;
5718        break;
5719        }
5720      }
5721      if (canMove) {
5722        if (changeUp*djBasic>1.0e-12||fabs(changeUp)<1.0e-8) {
5723          // move
5724          for (i = 0;i<rowLength[iRow];i++) {
5725            int iColumn = thisIndices[i];
5726            double alpha = thisElements[i];
5727            dj[iColumn] += djBasic * alpha;
5728          }
5729          dj[iSequence]=0.0;
5730#ifdef CLP_DEBUG
5731          {
5732            double value5=0.0;
5733            int i;
5734            for (i=0;i<numberRows_+numberColumns_;i++) {
5735              if (dj[i]<-1.0e-6)
5736                value5 += dj[i]*upper_[i];
5737              else if (dj[i] >1.0e-6)
5738                value5 += dj[i]*lower_[i];
5739            }
5740            printf("Values objective Value after row %d old dj %g %g\n",
5741                   iRow,djBasic,value5);
5742          }
5743#endif
5744        }
5745      }
5746    }     
5747  }
5748  delete [] scaled;
5749}
5750int
5751ClpSimplexDual::nextSuperBasic()
5752{
5753  if (firstFree_>=0) {
5754    int returnValue=firstFree_;
5755    int iColumn=firstFree_+1;
5756    for (;iColumn<numberRows_+numberColumns_;iColumn++) {
5757      if (getStatus(iColumn)==isFree) 
5758        if (fabs(dj_[iColumn])>1.0e2*dualTolerance_) 
5759          break;
5760    }
5761    firstFree_ = iColumn;
5762    if (firstFree_==numberRows_+numberColumns_)
5763      firstFree_=-1;
5764    return returnValue;
5765  } else {
5766    return -1;
5767  }
5768}
5769void 
5770ClpSimplexDual::resetFakeBounds()
5771{
5772  double dummyChangeCost=0.0;
5773  changeBounds(true,rowArray_[2],dummyChangeCost);
5774  // throw away change
5775  for (int i=0;i<4;i++) 
5776    rowArray_[i]->clear();
5777}
Note: See TracBrowser for help on using the repository browser.