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

Last change on this file since 866 was 866, checked in by forrest, 14 years ago

fix perturbation

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