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

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

for integer sync

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 173.4 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4
5/* Notes on implementation of dual simplex algorithm.
6
7   When dual feasible:
8
9   If primal feasible, we are optimal.  Otherwise choose an infeasible
10   basic variable to leave basis (normally going to nearest bound) (B).  We
11   now need to find an incoming variable which will leave problem
12   dual feasible so we get the row of the tableau corresponding to
13   the basic variable (with the correct sign depending if basic variable
14   above or below feasibility region - as that affects whether reduced
15   cost on outgoing variable has to be positive or negative).
16
17   We now perform a ratio test to determine which incoming variable will
18   preserve dual feasibility (C).  If no variable found then problem
19   is infeasible (in primal sense).  If there is a variable, we then
20   perform pivot and repeat.  Trivial?
21
22   -------------------------------------------
23
24   A) How do we get dual feasible?  If all variables have bounds then
25   it is trivial to get feasible by putting non-basic variables to
26   correct bounds.  OSL did not have a phase 1/phase 2 approach but
27   instead effectively put fake bounds on variables and this is the
28   approach here, although I had hoped to make it cleaner.
29
30   If there is a weight of X on getting dual feasible:
31     Non-basic variables with negative reduced costs are put to
32     lesser of their upper bound and their lower bound + X.
33     Similarly, mutatis mutandis, for positive reduced costs.
34
35   Free variables should normally be in basis, otherwise I have
36   coding which may be able to come out (and may not be correct).
37
38   In OSL, this weight was changed heuristically, here at present
39   it is only increased if problem looks finished.  If problem is
40   feasible I check for unboundedness.  If not unbounded we
41   could play with going into primal.  As long as weights increase
42   any algorithm would be finite.
43   
44   B) Which outgoing variable to choose is a virtual base class.
45   For difficult problems steepest edge is preferred while for
46   very easy (large) problems we will need partial scan.
47
48   C) Sounds easy, but this is hardest part of algorithm.
49      1) Instead of stopping at first choice, we may be able
50      to flip that variable to other bound and if objective
51      still improving choose again.  These mini iterations can
52      increase speed by orders of magnitude but we may need to
53      go to more of a bucket choice of variable rather than looking
54      at them one by one (for speed).
55      2) Accuracy.  Reduced costs may be of wrong sign but less than
56      tolerance.  Pivoting on these makes objective go backwards.
57      OSL modified cost so a zero move was made, Gill et al
58      (in primal analogue) modified so a strictly positive move was
59      made.  It is not quite as neat in dual but that is what we
60      try and do.  The two problems are that re-factorizations can
61      change reduced costs above and below tolerances and that when
62      finished we need to reset costs and try again.
63      3) Degeneracy.  Gill et al helps but may not be enough.  We
64      may need more.  Also it can improve speed a lot if we perturb
65      the costs significantly. 
66
67  References:
68     Forrest and Goldfarb, Steepest-edge simplex algorithms for
69       linear programming - Mathematical Programming 1992
70     Forrest and Tomlin, Implementing the simplex method for
71       the Optimization Subroutine Library - IBM Systems Journal 1992
72     Gill, Murray, Saunders, Wright A Practical Anti-Cycling
73       Procedure for Linear and Nonlinear Programming SOL report 1988
74
75
76  TODO:
77 
78  a) Better recovery procedures.  At present I never check on forward
79     progress.  There is checkpoint/restart with reducing
80     re-factorization frequency, but this is only on singular
81     factorizations.
82  b) Fast methods for large easy problems (and also the option for
83     the code to automatically choose which method).
84  c) We need to be able to stop in various ways for OSI - this
85     is fairly easy.
86
87 */
88
89
90#include "CoinPragma.hpp"
91
92#include <math.h>
93
94#include "CoinHelperFunctions.hpp"
95#include "ClpHelperFunctions.hpp"
96#include "ClpSimplexDual.hpp"
97#include "ClpEventHandler.hpp"
98#include "ClpFactorization.hpp"
99#include "CoinPackedMatrix.hpp"
100#include "CoinIndexedVector.hpp"
101#include "ClpDualRowDantzig.hpp"
102#include "ClpMessage.hpp"
103#include <cfloat>
104#include <cassert>
105#include <string>
106#include <stdio.h>
107#include <iostream>
108//#define CLP_DEBUG 1
109// To force to follow another run put logfile name here and define
110//#define FORCE_FOLLOW
111#ifdef FORCE_FOLLOW
112static FILE * fpFollow=NULL;
113static char * forceFile="old.log";
114static int force_in=-1;
115static int force_out=-1;
116static int force_iteration=0;
117#endif
118//#define VUB
119#ifdef VUB
120extern int * vub;
121extern int * toVub;
122extern int * nextDescendent;
123#endif
124// dual
125
126  /* *** Method
127     This is a vanilla version of dual simplex.
128
129     It tries to be a single phase approach with a weight of 1.0 being
130     given to getting optimal and a weight of dualBound_ being
131     given to getting dual feasible.  In this version I have used the
132     idea that this weight can be thought of as a fake bound.  If the
133     distance between the lower and upper bounds on a variable is less
134     than the feasibility weight then we are always better off flipping
135     to other bound to make dual feasible.  If the distance is greater
136     then we make up a fake bound dualBound_ away from one bound.
137     If we end up optimal or primal infeasible, we check to see if
138     bounds okay.  If so we have finished, if not we increase dualBound_
139     and continue (after checking if unbounded). I am undecided about
140     free variables - there is coding but I am not sure about it.  At
141     present I put them in basis anyway.
142
143     The code is designed to take advantage of sparsity so arrays are
144     seldom zeroed out from scratch or gone over in their entirety.
145     The only exception is a full scan to find outgoing variable.  This
146     will be changed to keep an updated list of infeasibilities (or squares
147     if steepest edge).  Also on easy problems we don't need full scan - just
148     pick first reasonable.
149
150     One problem is how to tackle degeneracy and accuracy.  At present
151     I am using the modification of costs which I put in OSL and which was
152     extended by Gill et al.  I am still not sure of the exact details.
153
154     The flow of dual is three while loops as follows:
155
156     while (not finished) {
157
158       while (not clean solution) {
159
160          Factorize and/or clean up solution by flipping variables so
161          dual feasible.  If looks finished check fake dual bounds.
162          Repeat until status is iterating (-1) or finished (0,1,2)
163
164       }
165
166       while (status==-1) {
167
168         Iterate until no pivot in or out or time to re-factorize.
169
170         Flow is:
171
172         choose pivot row (outgoing variable).  if none then
173         we are primal feasible so looks as if done but we need to
174         break and check bounds etc.
175
176         Get pivot row in tableau
177
178         Choose incoming column.  If we don't find one then we look
179         primal infeasible so break and check bounds etc.  (Also the
180         pivot tolerance is larger after any iterations so that may be
181         reason)
182
183         If we do find incoming column, we may have to adjust costs to
184         keep going forwards (anti-degeneracy).  Check pivot will be stable
185         and if unstable throw away iteration (we will need to implement
186         flagging of basic variables sometime) and break to re-factorize.
187         If minor error re-factorize after iteration.
188
189         Update everything (this may involve flipping variables to stay
190         dual feasible.
191
192       }
193
194     }
195
196     At present we never check we are going forwards.  I overdid that in
197     OSL so will try and make a last resort.
198
199     Needs partial scan pivot out option.
200     Needs dantzig, uninitialized and full steepest edge options (can still
201     use partial scan)
202
203     May need other anti-degeneracy measures, especially if we try and use
204     loose tolerances as a way to solve in fewer iterations.
205
206     I like idea of dynamic scaling.  This gives opportunity to decouple
207     different implications of scaling for accuracy, iteration count and
208     feasibility tolerance.
209
210  */
211#define CLEAN_FIXED 0
212// Startup part of dual (may be extended to other algorithms)
213int 
214ClpSimplexDual::startupSolve(int ifValuesPass,double * saveDuals,int startFinishOptions)
215{
216  // If values pass then save given duals round check solution
217  // sanity check
218  // initialize - no values pass and algorithm_ is -1
219  // put in standard form (and make row copy)
220  // create modifiable copies of model rim and do optional scaling
221  // If problem looks okay
222  // Do initial factorization
223  // If user asked for perturbation - do it
224  if (!startup(0,startFinishOptions)) {
225    // looks okay
226    // Superbasic variables not allowed
227    // If values pass then scale pi
228    if (ifValuesPass) {
229      if (problemStatus_&&perturbation_<100) 
230        perturb();
231      int i;
232      if (scalingFlag_>0) {
233        for (i=0;i<numberRows_;i++) {
234          dual_[i] = saveDuals[i]/rowScale_[i];
235        }
236      } else {
237        CoinMemcpyN(saveDuals,numberRows_,dual_);
238      }
239      // now create my duals
240      for (i=0;i<numberRows_;i++) {
241        // slack
242        double value = dual_[i];
243        value += rowObjectiveWork_[i];
244        saveDuals[i+numberColumns_]=value;
245      }
246      CoinMemcpyN(objectiveWork_,numberColumns_,saveDuals);
247      transposeTimes(-1.0,dual_,saveDuals);
248      // make reduced costs okay
249      for (i=0;i<numberColumns_;i++) {
250        if (getStatus(i)==atLowerBound) {
251          if (saveDuals[i]<0.0) {
252            //if (saveDuals[i]<-1.0e-3)
253            //printf("bad dj at lb %d %g\n",i,saveDuals[i]);
254            saveDuals[i]=0.0;
255          }
256        } else if (getStatus(i)==atUpperBound) {
257          if (saveDuals[i]>0.0) {
258            //if (saveDuals[i]>1.0e-3)
259            //printf("bad dj at ub %d %g\n",i,saveDuals[i]);
260            saveDuals[i]=0.0;
261          }
262        }
263      }
264      CoinMemcpyN(saveDuals,(numberColumns_+numberRows_),dj_);
265      // set up possible ones
266      for (i=0;i<numberRows_+numberColumns_;i++)
267        clearPivoted(i);
268      int iRow;
269      for (iRow=0;iRow<numberRows_;iRow++) {
270        int iPivot=pivotVariable_[iRow];
271        if (fabs(saveDuals[iPivot])>dualTolerance_) {
272          if (getStatus(iPivot)!=isFree) 
273            setPivoted(iPivot);
274        }
275      }
276    } else if ((specialOptions_&1024)!=0&&CLEAN_FIXED) {
277      // set up possible ones
278      for (int i=0;i<numberRows_+numberColumns_;i++)
279        clearPivoted(i);
280      int iRow;
281      for (iRow=0;iRow<numberRows_;iRow++) {
282        int iPivot=pivotVariable_[iRow];
283        if (iPivot<numberColumns_&&lower_[iPivot]==upper_[iPivot]) {
284          setPivoted(iPivot);
285        }
286      }
287    } 
288
289    double objectiveChange;
290    numberFake_ =0; // Number of variables at fake bounds
291    numberChanged_ =0; // Number of variables with changed costs
292    changeBounds(true,NULL,objectiveChange);
293   
294    if (!ifValuesPass) {
295      // Check optimal
296      if (!numberDualInfeasibilities_&&!numberPrimalInfeasibilities_)
297        problemStatus_=0;
298    }
299    if (problemStatus_<0&&perturbation_<100) {
300      bool inCbcOrOther = (specialOptions_&0x03000000)!=0;
301      if (!inCbcOrOther)
302        perturb();
303      // Can't get here if values pass
304      gutsOfSolution(NULL,NULL);
305      if (handler_->logLevel()>2) {
306        handler_->message(CLP_SIMPLEX_STATUS,messages_)
307          <<numberIterations_<<objectiveValue();
308        handler_->printing(sumPrimalInfeasibilities_>0.0)
309          <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
310        handler_->printing(sumDualInfeasibilities_>0.0)
311          <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
312        handler_->printing(numberDualInfeasibilitiesWithoutFree_
313                           <numberDualInfeasibilities_)
314                             <<numberDualInfeasibilitiesWithoutFree_;
315        handler_->message()<<CoinMessageEol;
316      }
317      if (inCbcOrOther) {
318        if (numberPrimalInfeasibilities_) {
319          perturb();
320          if (perturbation_>=101)
321            computeDuals(NULL);
322          //gutsOfSolution(NULL,NULL);
323        } else if (numberDualInfeasibilities_) {
324          problemStatus_=10;
325          return 1; // to primal
326        }
327      }
328    }
329    return 0;
330  } else {
331    return 1;
332  }
333}
334void 
335ClpSimplexDual::finishSolve(int startFinishOptions)
336{
337  assert (problemStatus_||!sumPrimalInfeasibilities_);
338
339  // clean up
340  finish(startFinishOptions);
341}
342void 
343ClpSimplexDual::gutsOfDual(int ifValuesPass,double * & saveDuals,int initialStatus,
344                           ClpDataSave & data)
345{
346  int lastCleaned=0; // last time objective or bounds cleaned up
347 
348  // This says whether to restore things etc
349  // startup will have factorized so can skip
350  int factorType=0;
351  // Start check for cycles
352  progress_->startCheck();
353  // Say change made on first iteration
354  changeMade_=1;
355  /*
356    Status of problem:
357    0 - optimal
358    1 - infeasible
359    2 - unbounded
360    -1 - iterating
361    -2 - factorization wanted
362    -3 - redo checking without factorization
363    -4 - looks infeasible
364  */
365  while (problemStatus_<0) {
366    int iRow, iColumn;
367    // clear
368    for (iRow=0;iRow<4;iRow++) {
369      rowArray_[iRow]->clear();
370    }   
371   
372    for (iColumn=0;iColumn<2;iColumn++) {
373      columnArray_[iColumn]->clear();
374    }   
375   
376    // give matrix (and model costs and bounds a chance to be
377    // refreshed (normally null)
378    matrix_->refresh(this);
379    // If getting nowhere - why not give it a kick
380    // does not seem to work too well - do some more work
381    if (perturbation_<101&&numberIterations_>2*(numberRows_+numberColumns_)
382        &&initialStatus!=10) {
383      perturb();
384      // Can't get here if values pass
385      gutsOfSolution(NULL,NULL);
386      if (handler_->logLevel()>2) {
387        handler_->message(CLP_SIMPLEX_STATUS,messages_)
388          <<numberIterations_<<objectiveValue();
389        handler_->printing(sumPrimalInfeasibilities_>0.0)
390          <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
391        handler_->printing(sumDualInfeasibilities_>0.0)
392          <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
393        handler_->printing(numberDualInfeasibilitiesWithoutFree_
394                           <numberDualInfeasibilities_)
395                             <<numberDualInfeasibilitiesWithoutFree_;
396        handler_->message()<<CoinMessageEol;
397      }
398    }
399    // see if in Cbc etc
400    bool inCbcOrOther = (specialOptions_&0x03000000)!=0;
401#if 0
402    bool gotoPrimal=false;
403    if (inCbcOrOther&&numberIterations_>disasterArea_+numberRows_&&
404        numberDualInfeasibilitiesWithoutFree_&&largestDualError_>1.0e-1) {
405      if (!disasterArea_) {
406        printf("trying all slack\n");
407        // try all slack basis
408        allSlackBasis(true);
409        disasterArea_=2*numberRows_;
410      } else {
411        printf("going to primal\n");
412        // go to primal
413        gotoPrimal=true;
414        allSlackBasis(true);
415      }
416    }
417#endif
418    bool disaster=false;
419    if (disasterArea_&&inCbcOrOther&&disasterArea_->check()) {
420      disasterArea_->saveInfo();
421      disaster=true;
422    }
423    // may factorize, checks if problem finished
424    statusOfProblemInDual(lastCleaned,factorType,saveDuals,data,
425                          ifValuesPass);
426    if (disaster)
427      problemStatus_=3;
428    // If values pass then do easy ones on first time
429    if (ifValuesPass&&
430        progress_->lastIterationNumber(0)<0&&saveDuals) {
431      doEasyOnesInValuesPass(saveDuals);
432    }
433   
434    // Say good factorization
435    factorType=1;
436    if (data.sparseThreshold_) {
437      // use default at present
438      factorization_->sparseThreshold(0);
439      factorization_->goSparse();
440    }
441   
442    // exit if victory declared
443    if (problemStatus_>=0)
444      break;
445   
446    // test for maximum iterations
447    if (hitMaximumIterations()||(ifValuesPass==2&&!saveDuals)) {
448      problemStatus_=3;
449      break;
450    }
451    if (ifValuesPass&&!saveDuals) {
452      // end of values pass
453      ifValuesPass=0;
454      int status = eventHandler_->event(ClpEventHandler::endOfValuesPass);
455      if (status>=0) {
456        problemStatus_=5;
457        secondaryStatus_=ClpEventHandler::endOfValuesPass;
458        break;
459      }
460    }
461    // Check event
462    {
463      int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
464      if (status>=0) {
465        problemStatus_=5;
466        secondaryStatus_=ClpEventHandler::endOfFactorization;
467        break;
468      }
469    }
470      // Do iterations
471    whileIterating(saveDuals,ifValuesPass);
472  }
473}
474int 
475ClpSimplexDual::dual(int ifValuesPass,int startFinishOptions)
476{
477  algorithm_ = -1;
478 
479  // save data
480  ClpDataSave data = saveData();
481  double * saveDuals = NULL;
482  if (ifValuesPass) {
483    saveDuals = new double [numberRows_+numberColumns_];
484    CoinMemcpyN(dual_,numberRows_,saveDuals);
485  }
486  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        if (printOut)
4607          printf("col %d cost %g change %g\n",iColumn,objectiveWork_[iColumn],value);
4608        objectiveWork_[iColumn] += value;
4609      }
4610    }
4611  }
4612  handler_->message(CLP_SIMPLEX_PERTURB,messages_)
4613    <<100.0*maximumFraction<<perturbation<<largest<<100.0*largestPerCent<<largestZero
4614    <<CoinMessageEol;
4615  // and zero changes
4616  //int nTotal = numberRows_+numberColumns_;
4617  //CoinZeroN(cost_+nTotal,nTotal);
4618  // say perturbed
4619  perturbation_=101;
4620
4621}
4622/* For strong branching.  On input lower and upper are new bounds
4623   while on output they are change in objective function values
4624   (>1.0e50 infeasible).
4625   Return code is 0 if nothing interesting, -1 if infeasible both
4626   ways and +1 if infeasible one way (check values to see which one(s))
4627   Returns -2 if bad factorization
4628*/
4629int ClpSimplexDual::strongBranching(int numberVariables,const int * variables,
4630                                    double * newLower, double * newUpper,
4631                                    double ** outputSolution,
4632                                    int * outputStatus, int * outputIterations,
4633                                    bool stopOnFirstInfeasible,
4634                                    bool alwaysFinish,
4635                                    int startFinishOptions)
4636{
4637  int i;
4638  int returnCode=0;
4639  double saveObjectiveValue = objectiveValue_;
4640  algorithm_ = -1;
4641
4642  //scaling(false);
4643
4644  // put in standard form (and make row copy)
4645  // create modifiable copies of model rim and do optional scaling
4646  createRim(7+8+16+32,true,startFinishOptions);
4647
4648  // change newLower and newUpper if scaled
4649
4650  // Do initial factorization
4651  // and set certain stuff
4652  // We can either set increasing rows so ...IsBasic gives pivot row
4653  // or we can just increment iBasic one by one
4654  // for now let ...iBasic give pivot row
4655  int useFactorization=false;
4656  if ((startFinishOptions&2)!=0&&(whatsChanged_&(2+512))==2+512)
4657    useFactorization=true; // Keep factorization if possible
4658  // switch off factorization if bad
4659  if (pivotVariable_[0]<0)
4660    useFactorization=false;
4661  if (!useFactorization||factorization_->numberRows()!=numberRows_) {
4662    useFactorization = false;
4663    factorization_->increasingRows(2);
4664    // row activities have negative sign
4665    factorization_->slackValue(-1.0);
4666    factorization_->zeroTolerance(1.0e-13);
4667
4668    int factorizationStatus = internalFactorize(0);
4669    if (factorizationStatus<0) {
4670      // some error
4671      // we should either debug or ignore
4672#ifndef NDEBUG
4673      printf("***** ClpDual strong branching factorization error - debug\n");
4674#endif
4675      return -2;
4676    } else if (factorizationStatus&&factorizationStatus<=numberRows_) {
4677      handler_->message(CLP_SINGULARITIES,messages_)
4678        <<factorizationStatus
4679        <<CoinMessageEol;
4680    }
4681  }
4682  // save stuff
4683  ClpFactorization saveFactorization(*factorization_);
4684  // save basis and solution
4685  double * saveSolution = new double[numberRows_+numberColumns_];
4686  CoinMemcpyN(solution_,
4687         numberRows_+numberColumns_,saveSolution);
4688  unsigned char * saveStatus = 
4689    new unsigned char [numberRows_+numberColumns_];
4690  CoinMemcpyN(status_,numberColumns_+numberRows_,saveStatus);
4691  // save bounds as createRim makes clean copies
4692  double * saveLower = new double[numberRows_+numberColumns_];
4693  CoinMemcpyN(lower_,
4694         numberRows_+numberColumns_,saveLower);
4695  double * saveUpper = new double[numberRows_+numberColumns_];
4696  CoinMemcpyN(upper_,
4697         numberRows_+numberColumns_,saveUpper);
4698  double * saveObjective = new double[numberRows_+numberColumns_];
4699  CoinMemcpyN(cost_,
4700         numberRows_+numberColumns_,saveObjective);
4701  int * savePivot = new int [numberRows_];
4702  CoinMemcpyN(pivotVariable_, numberRows_,savePivot);
4703  // need to save/restore weights.
4704
4705  int iSolution = 0;
4706  for (i=0;i<numberVariables;i++) {
4707    int iColumn = variables[i];
4708    double objectiveChange;
4709    double saveBound;
4710   
4711    // try down
4712   
4713    saveBound = columnUpper_[iColumn];
4714    // external view - in case really getting optimal
4715    columnUpper_[iColumn] = newUpper[i];
4716    if (scalingFlag_<=0) 
4717      upper_[iColumn] = newUpper[i]*rhsScale_;
4718    else 
4719      upper_[iColumn] = (newUpper[i]/columnScale_[iColumn])*rhsScale_; // scale
4720    // Start of fast iterations
4721    int status = fastDual(alwaysFinish);
4722    CoinAssert (problemStatus_||objectiveValue_<1.0e50);
4723    // make sure plausible
4724    double obj = CoinMax(objectiveValue_,saveObjectiveValue);
4725    if (status&&problemStatus_!=3) {
4726      // not finished - might be optimal
4727      checkPrimalSolution(rowActivityWork_,columnActivityWork_);
4728      double limit = 0.0;
4729      getDblParam(ClpDualObjectiveLimit, limit);
4730      if (!numberPrimalInfeasibilities_&&obj<limit) { 
4731        problemStatus_=0;
4732      } 
4733      status=problemStatus_;
4734    }
4735    if (status||(problemStatus_==0&&!isDualObjectiveLimitReached())) {
4736      objectiveChange = obj-saveObjectiveValue;
4737    } else {
4738      objectiveChange = 1.0e100;
4739      status=1;
4740    }
4741    if (problemStatus_==3)
4742      status=2;
4743
4744    if (scalingFlag_<=0) {
4745      CoinMemcpyN(solution_,numberColumns_,outputSolution[iSolution]);
4746    } else {
4747      int j;
4748      double * sol = outputSolution[iSolution];
4749      for (j=0;j<numberColumns_;j++) 
4750        sol[j] = solution_[j]*columnScale_[j];
4751    }
4752    outputStatus[iSolution]=status;
4753    outputIterations[iSolution]=numberIterations_;
4754    iSolution++;
4755    // restore
4756    CoinMemcpyN(saveSolution,
4757            numberRows_+numberColumns_,solution_);
4758    CoinMemcpyN(saveStatus,numberColumns_+numberRows_,status_);
4759    CoinMemcpyN(saveLower,
4760            numberRows_+numberColumns_,lower_);
4761    CoinMemcpyN(saveUpper,
4762            numberRows_+numberColumns_,upper_);
4763    CoinMemcpyN(saveObjective,
4764            numberRows_+numberColumns_,cost_);
4765    columnUpper_[iColumn] = saveBound;
4766    CoinMemcpyN(savePivot, numberRows_,pivotVariable_);
4767    delete factorization_;
4768    factorization_ = new ClpFactorization(saveFactorization);
4769
4770    newUpper[i]=objectiveChange;
4771#ifdef CLP_DEBUG
4772    printf("down on %d costs %g\n",iColumn,objectiveChange);
4773#endif
4774         
4775    // try up
4776   
4777    saveBound = columnLower_[iColumn];
4778    // external view - in case really getting optimal
4779    columnLower_[iColumn] = newLower[i];
4780    if (scalingFlag_<=0) 
4781      lower_[iColumn] = newLower[i]*rhsScale_;
4782    else 
4783      lower_[iColumn] = (newLower[i]/columnScale_[iColumn])*rhsScale_; // scale
4784    // Start of fast iterations
4785    status = fastDual(alwaysFinish);
4786    // make sure plausible
4787    obj = CoinMax(objectiveValue_,saveObjectiveValue);
4788    if (status&&problemStatus_!=3) {
4789      // not finished - might be optimal
4790      checkPrimalSolution(rowActivityWork_,columnActivityWork_);
4791      double limit = 0.0;
4792      getDblParam(ClpDualObjectiveLimit, limit);
4793      if (!numberPrimalInfeasibilities_&&obj< limit) { 
4794        problemStatus_=0;
4795      } 
4796      status=problemStatus_;
4797    }
4798    if (status||(problemStatus_==0&&!isDualObjectiveLimitReached())) {
4799      objectiveChange = obj-saveObjectiveValue;
4800    } else {
4801      objectiveChange = 1.0e100;
4802      status=1;
4803    }
4804    if (problemStatus_==3)
4805      status=2;
4806    if (scalingFlag_<=0) {
4807      CoinMemcpyN(solution_,numberColumns_,outputSolution[iSolution]);
4808    } else {
4809      int j;
4810      double * sol = outputSolution[iSolution];
4811      for (j=0;j<numberColumns_;j++) 
4812        sol[j] = solution_[j]*columnScale_[j];
4813    }
4814    outputStatus[iSolution]=status;
4815    outputIterations[iSolution]=numberIterations_;
4816    iSolution++;
4817
4818    // restore
4819    CoinMemcpyN(saveSolution,
4820            numberRows_+numberColumns_,solution_);
4821    CoinMemcpyN(saveStatus,numberColumns_+numberRows_,status_);
4822    CoinMemcpyN(saveLower,
4823            numberRows_+numberColumns_,lower_);
4824    CoinMemcpyN(saveUpper,
4825            numberRows_+numberColumns_,upper_);
4826    CoinMemcpyN(saveObjective,
4827            numberRows_+numberColumns_,cost_);
4828    columnLower_[iColumn] = saveBound;
4829    CoinMemcpyN(savePivot, numberRows_,pivotVariable_);
4830    delete factorization_;
4831    factorization_ = new ClpFactorization(saveFactorization);
4832
4833    newLower[i]=objectiveChange;
4834#ifdef CLP_DEBUG
4835    printf("up on %d costs %g\n",iColumn,objectiveChange);
4836#endif
4837         
4838    /* Possibilities are:
4839       Both sides feasible - store
4840       Neither side feasible - set objective high and exit if desired
4841       One side feasible - change bounds and resolve
4842    */
4843    if (newUpper[i]<1.0e100) {
4844      if(newLower[i]<1.0e100) {
4845        // feasible - no action
4846      } else {
4847        // up feasible, down infeasible
4848        returnCode=1;
4849        if (stopOnFirstInfeasible)
4850          break;
4851      }
4852    } else {
4853      if(newLower[i]<1.0e100) {
4854        // down feasible, up infeasible
4855        returnCode=1;
4856        if (stopOnFirstInfeasible)
4857          break;
4858      } else {
4859        // neither side feasible
4860        returnCode=-1;
4861        break;
4862      }
4863    }
4864  }
4865  delete [] saveSolution;
4866  delete [] saveLower;
4867  delete [] saveUpper;
4868  delete [] saveObjective;
4869  delete [] saveStatus;
4870  delete [] savePivot;
4871  if ((startFinishOptions&1)==0) {
4872    deleteRim(1);
4873    whatsChanged_=0;
4874  } else {
4875    // Original factorization will have been put back by last loop
4876    //delete factorization_;
4877    //factorization_ = new ClpFactorization(saveFactorization);
4878    deleteRim(0);
4879    // mark all as current
4880    whatsChanged_ = 0xffff;
4881  }
4882  objectiveValue_ = saveObjectiveValue;
4883  return returnCode;
4884}
4885// treat no pivot as finished (unless interesting)
4886int ClpSimplexDual::fastDual(bool alwaysFinish)
4887{
4888  algorithm_ = -1;
4889  secondaryStatus_=0;
4890  // Say in fast dual
4891  specialOptions_ |= 16384;
4892  //handler_->setLogLevel(63);
4893  // save data
4894  ClpDataSave data = saveData();
4895  dualTolerance_=dblParam_[ClpDualTolerance];
4896  primalTolerance_=dblParam_[ClpPrimalTolerance];
4897
4898  // save dual bound
4899  double saveDualBound = dualBound_;
4900
4901  double objectiveChange;
4902  // for dual we will change bounds using dualBound_
4903  // for this we need clean basis so it is after factorize
4904  gutsOfSolution(NULL,NULL);
4905  numberFake_ =0; // Number of variables at fake bounds
4906  numberChanged_ =0; // Number of variables with changed costs
4907  changeBounds(true,NULL,objectiveChange);
4908
4909  problemStatus_ = -1;
4910  numberIterations_=0;
4911  factorization_->sparseThreshold(0);
4912  factorization_->goSparse();
4913
4914  int lastCleaned=0; // last time objective or bounds cleaned up
4915
4916  // number of times we have declared optimality
4917  numberTimesOptimal_=0;
4918
4919  // This says whether to restore things etc
4920  int factorType=0;
4921  /*
4922    Status of problem:
4923    0 - optimal
4924    1 - infeasible
4925    2 - unbounded
4926    -1 - iterating
4927    -2 - factorization wanted
4928    -3 - redo checking without factorization
4929    -4 - looks infeasible
4930
4931    BUT also from whileIterating return code is:
4932
4933   -1 iterations etc
4934   -2 inaccuracy
4935   -3 slight inaccuracy (and done iterations)
4936   +0 looks optimal (might be unbounded - but we will investigate)
4937   +1 looks infeasible
4938   +3 max iterations
4939
4940  */
4941
4942  int returnCode = 0;
4943
4944  int iRow,iColumn;
4945  while (problemStatus_<0) {
4946    // clear
4947    for (iRow=0;iRow<4;iRow++) {
4948      rowArray_[iRow]->clear();
4949    }   
4950   
4951    for (iColumn=0;iColumn<2;iColumn++) {
4952      columnArray_[iColumn]->clear();
4953    }   
4954
4955    // give matrix (and model costs and bounds a chance to be
4956    // refreshed (normally null)
4957    matrix_->refresh(this);
4958    // may factorize, checks if problem finished
4959    // should be able to speed this up on first time
4960    statusOfProblemInDual(lastCleaned,factorType,NULL,data,0);
4961
4962    // Say good factorization
4963    factorType=1;
4964
4965    // Do iterations
4966    if (problemStatus_<0) {
4967      double * givenPi=NULL;
4968      returnCode = whileIterating(givenPi,0);
4969      if ((!alwaysFinish&&returnCode<0)||returnCode==3) {
4970        if (returnCode!=3)
4971          assert (problemStatus_<0);
4972        returnCode=1;
4973        problemStatus_ = 3; 
4974        // can't say anything interesting - might as well return
4975#ifdef CLP_DEBUG
4976        printf("returning from fastDual after %d iterations with code %d\n",
4977               numberIterations_,returnCode);
4978#endif
4979        break;
4980      }
4981      returnCode=0;
4982    }
4983  }
4984
4985  // clear
4986  for (iRow=0;iRow<4;iRow++) {
4987    rowArray_[iRow]->clear();
4988  }   
4989 
4990  for (iColumn=0;iColumn<2;iColumn++) {
4991    columnArray_[iColumn]->clear();
4992  }   
4993  // Say not in fast dual
4994  specialOptions_ &= ~16384;
4995  assert(!numberFake_||((specialOptions_&(2048|4096))!=0&&dualBound_>1.0e8)
4996         ||returnCode||problemStatus_); // all bounds should be okay
4997  // Restore any saved stuff
4998  restoreData(data);
4999  dualBound_ = saveDualBound;
5000  return returnCode;
5001}
5002// This does first part of StrongBranching
5003ClpFactorization * 
5004ClpSimplexDual::setupForStrongBranching(char * arrays, int numberRows, int numberColumns)
5005{
5006  algorithm_ = -1;
5007  // put in standard form (and make row copy)
5008  // create modifiable copies of model rim and do optional scaling
5009  int startFinishOptions;
5010  /*  COIN_CLP_VETTED
5011      Looks safe for Cbc
5012  */
5013  bool safeOption = ((specialOptions_&COIN_CBC_USING_CLP)!=0);
5014  if((specialOptions_&4096)==0||(auxiliaryModel_&&!safeOption)) {
5015    startFinishOptions=0;
5016  } else {
5017    startFinishOptions=1+2+4;
5018  }
5019  createRim(7+8+16+32,true,startFinishOptions);
5020  // Do initial factorization
5021  // and set certain stuff
5022  // We can either set increasing rows so ...IsBasic gives pivot row
5023  // or we can just increment iBasic one by one
5024  // for now let ...iBasic give pivot row
5025  bool useFactorization=false;
5026  if ((startFinishOptions&2)!=0&&(whatsChanged_&(2+512))==2+512) {
5027    useFactorization=true; // Keep factorization if possible
5028    // switch off factorization if bad
5029    if (pivotVariable_[0]<0||factorization_->numberRows()!=numberRows_) 
5030      useFactorization=false;
5031  }
5032  if (!useFactorization) {
5033    factorization_->increasingRows(2);
5034    // row activities have negative sign
5035    factorization_->slackValue(-1.0);
5036    factorization_->zeroTolerance(1.0e-13);
5037
5038    int factorizationStatus = internalFactorize(0);
5039    if (factorizationStatus<0) {
5040      // some error
5041      // we should either debug or ignore
5042#ifndef NDEBUG
5043      printf("***** ClpDual strong branching factorization error - debug\n");
5044#endif
5045    } else if (factorizationStatus&&factorizationStatus<=numberRows_) {
5046      handler_->message(CLP_SINGULARITIES,messages_)
5047        <<factorizationStatus
5048        <<CoinMessageEol;
5049    }
5050  }
5051  double * arrayD = (double *) arrays;
5052  arrayD[0]=objectiveValue()*optimizationDirection_;
5053  double * saveSolution = arrayD+1;
5054  double * saveLower = saveSolution + (numberRows+numberColumns);
5055  double * saveUpper = saveLower + (numberRows+numberColumns);
5056  double * saveObjective = saveUpper + (numberRows+numberColumns);
5057  double * saveLowerOriginal = saveObjective + (numberRows+numberColumns);
5058  double * saveUpperOriginal = saveLowerOriginal + numberColumns;
5059  arrayD = saveUpperOriginal + numberColumns;
5060  int * savePivot = (int *) arrayD;
5061  int * whichRow = savePivot+numberRows;
5062  int * whichColumn = whichRow + 3*numberRows;
5063  int * arrayI = whichColumn + 2*numberColumns;
5064  unsigned char * saveStatus = (unsigned char *) (arrayI+1);
5065  // save stuff
5066  // save basis and solution
5067  CoinMemcpyN(solution_,
5068          numberRows_+numberColumns_,saveSolution);
5069  CoinMemcpyN(status_,numberColumns_+numberRows_,saveStatus);
5070  CoinMemcpyN(lower_,
5071          numberRows_+numberColumns_,saveLower);
5072  CoinMemcpyN(upper_,
5073          numberRows_+numberColumns_,saveUpper);
5074  CoinMemcpyN(cost_,
5075          numberRows_+numberColumns_,saveObjective);
5076  CoinMemcpyN(pivotVariable_, numberRows_,savePivot);
5077  return new ClpFactorization(*factorization_);
5078}
5079// This cleans up after strong branching
5080void 
5081ClpSimplexDual::cleanupAfterStrongBranching()
5082{
5083  int startFinishOptions;
5084  /*  COIN_CLP_VETTED
5085      Looks safe for Cbc
5086  */
5087  bool safeOption = ((specialOptions_&COIN_CBC_USING_CLP)!=0);
5088  if((specialOptions_&4096)==0||(auxiliaryModel_&&!safeOption)) {
5089    startFinishOptions=0;
5090  } else {
5091    startFinishOptions=1+2+4;
5092  }
5093  if ((startFinishOptions&1)==0) {
5094    deleteRim(1);
5095    whatsChanged_=0;
5096  } else {
5097    // Original factorization will have been put back by last loop
5098    //delete factorization_;
5099    //factorization_ = new ClpFactorization(saveFactorization);
5100    deleteRim(0);
5101    // mark all as current
5102    whatsChanged_ = 0xffff;
5103  }
5104}
5105/* Checks number of variables at fake bounds.  This is used by fastDual
5106   so can exit gracefully before end */
5107int 
5108ClpSimplexDual::numberAtFakeBound()
5109{
5110  int iSequence;
5111  int numberFake=0;
5112 
5113  for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
5114    FakeBound bound = getFakeBound(iSequence);
5115    switch(getStatus(iSequence)) {
5116
5117    case basic:
5118      break;
5119    case isFree:
5120    case superBasic:
5121    case ClpSimplex::isFixed:
5122      assert (bound==noFake);
5123      break;
5124    case atUpperBound:
5125      if (bound==upperFake||bound==bothFake)
5126        numberFake++;
5127      break;
5128    case atLowerBound:
5129      if (bound==lowerFake||bound==bothFake)
5130        numberFake++;
5131      break;
5132    }
5133  }
5134  numberFake_ = numberFake;
5135  return numberFake;
5136}
5137/* Pivot out a variable and choose an incoing one.  Assumes dual
5138   feasible - will not go through a reduced cost. 
5139   Returns step length in theta
5140   Returns ray in ray_ (or NULL if no pivot)
5141   Return codes as before but -1 means no acceptable pivot
5142*/
5143int 
5144ClpSimplexDual::pivotResult()
5145{
5146  abort();
5147  return 0;
5148}
5149/*
5150   Row array has row part of pivot row
5151   Column array has column part.
5152   This is used in dual values pass
5153*/
5154void
5155ClpSimplexDual::checkPossibleValuesMove(CoinIndexedVector * rowArray,
5156                                        CoinIndexedVector * columnArray,
5157                                        double acceptablePivot)
5158{
5159  double * work;
5160  int number;
5161  int * which;
5162  int iSection;
5163
5164  double tolerance = dualTolerance_*1.001;
5165
5166  double thetaDown = 1.0e31;
5167  double changeDown ;
5168  double thetaUp = 1.0e31;
5169  double bestAlphaDown = acceptablePivot*0.99999;
5170  double bestAlphaUp = acceptablePivot*0.99999;
5171  int sequenceDown =-1;
5172  int sequenceUp = sequenceOut_;
5173
5174  double djBasic = dj_[sequenceOut_];
5175  if (djBasic>0.0) {
5176    // basic at lower bound so directionOut_ 1 and -1 in pivot row
5177    // dj will go to zero on other way
5178    thetaUp = djBasic;
5179    changeDown = -lower_[sequenceOut_];
5180  } else {
5181    // basic at upper bound so directionOut_ -1 and 1 in pivot row
5182    // dj will go to zero on other way
5183    thetaUp = -djBasic;
5184    changeDown = upper_[sequenceOut_];
5185  }
5186  bestAlphaUp = 1.0;
5187  int addSequence;
5188
5189  double alphaUp=0.0;
5190  double alphaDown=0.0;
5191
5192  for (iSection=0;iSection<2;iSection++) {
5193
5194    int i;
5195    if (!iSection) {
5196      work = rowArray->denseVector();
5197      number = rowArray->getNumElements();
5198      which = rowArray->getIndices();
5199      addSequence = numberColumns_;
5200    } else {
5201      work = columnArray->denseVector();
5202      number = columnArray->getNumElements();
5203      which = columnArray->getIndices();
5204      addSequence = 0;
5205    }
5206   
5207    for (i=0;i<number;i++) {
5208      int iSequence = which[i];
5209      int iSequence2 = iSequence + addSequence;
5210      double alpha;
5211      double oldValue;
5212      double value;
5213
5214      switch(getStatus(iSequence2)) {
5215         
5216      case basic: 
5217        break;
5218      case ClpSimplex::isFixed:
5219        alpha = work[i];
5220        changeDown += alpha*upper_[iSequence2];
5221        break;
5222      case isFree:
5223      case superBasic:
5224        alpha = work[i];
5225        // dj must be effectively zero as dual feasible
5226        if (fabs(alpha)>bestAlphaUp) {
5227          thetaDown = 0.0;
5228          thetaUp = 0.0;
5229          bestAlphaDown = fabs(alpha);
5230          bestAlphaUp = bestAlphaUp;
5231          sequenceDown =iSequence2;
5232          sequenceUp = sequenceDown;
5233          alphaUp = alpha;
5234          alphaDown = alpha;
5235        }
5236        break;
5237      case atUpperBound:
5238        alpha = work[i];
5239        oldValue = dj_[iSequence2];
5240        changeDown += alpha*upper_[iSequence2];
5241        if (alpha>=acceptablePivot) {
5242          // might do other way
5243          value = oldValue+thetaUp*alpha;
5244          if (value>-tolerance) {
5245            if (value>tolerance||fabs(alpha)>bestAlphaUp) {
5246              thetaUp = -oldValue/alpha;
5247              bestAlphaUp = fabs(alpha);
5248              sequenceUp = iSequence2;
5249              alphaUp = alpha;
5250            }
5251          }
5252        } else if (alpha<=-acceptablePivot) {
5253          // might do this way
5254          value = oldValue-thetaDown*alpha;
5255          if (value>-tolerance) {
5256            if (value>tolerance||fabs(alpha)>bestAlphaDown) {
5257              thetaDown = oldValue/alpha;
5258              bestAlphaDown = fabs(alpha);
5259              sequenceDown = iSequence2;
5260              alphaDown = alpha;
5261            }
5262          }
5263        }
5264        break;
5265      case atLowerBound:
5266        alpha = work[i];
5267        oldValue = dj_[iSequence2];
5268        changeDown += alpha*lower_[iSequence2];
5269        if (alpha<=-acceptablePivot) {
5270          // might do other way
5271          value = oldValue+thetaUp*alpha;
5272          if (value<tolerance) {
5273            if (value<-tolerance||fabs(alpha)>bestAlphaUp) {
5274              thetaUp = -oldValue/alpha;
5275              bestAlphaUp = fabs(alpha);
5276              sequenceUp = iSequence2;
5277              alphaUp = alpha;
5278            }
5279          }
5280        } else if (alpha>=acceptablePivot) {
5281          // might do this way
5282          value = oldValue-thetaDown*alpha;
5283          if (value<tolerance) {
5284            if (value<-tolerance||fabs(alpha)>bestAlphaDown) {
5285              thetaDown = oldValue/alpha;
5286              bestAlphaDown = fabs(alpha);
5287              sequenceDown = iSequence2;
5288              alphaDown = alpha;
5289            }
5290          }
5291        }
5292        break;
5293      }
5294    }
5295  }
5296  thetaUp *= -1.0;
5297  double changeUp = -thetaUp * changeDown;
5298  changeDown = -thetaDown*changeDown;
5299  if (CoinMax(fabs(thetaDown),fabs(thetaUp))<1.0e-8) {
5300    // largest
5301    if (fabs(alphaDown)<fabs(alphaUp)) {
5302      sequenceDown =-1;
5303    }
5304  }
5305  // choose
5306  sequenceIn_=-1;
5307  if (changeDown>changeUp&&sequenceDown>=0) {
5308    theta_ = thetaDown;
5309    if (fabs(changeDown)<1.0e30)
5310      sequenceIn_ = sequenceDown;
5311    alpha_ = alphaDown;
5312#ifdef CLP_DEBUG
5313    if ((handler_->logLevel()&32))
5314      printf("predicted way - dirout %d, change %g,%g theta %g\n",
5315             directionOut_,changeDown,changeUp,theta_);
5316#endif
5317  } else {
5318    theta_ = thetaUp;
5319    if (fabs(changeUp)<1.0e30)
5320      sequenceIn_ = sequenceUp;
5321    alpha_ = alphaUp;
5322    if (sequenceIn_!=sequenceOut_) {
5323#ifdef CLP_DEBUG
5324      if ((handler_->logLevel()&32))
5325        printf("opposite way - dirout %d, change %g,%g theta %g\n",
5326               directionOut_,changeDown,changeUp,theta_);
5327#endif
5328    } else {
5329#ifdef CLP_DEBUG
5330      if ((handler_->logLevel()&32))
5331        printf("opposite way to zero dj - dirout %d, change %g,%g theta %g\n",
5332               directionOut_,changeDown,changeUp,theta_);
5333#endif
5334    }
5335  }
5336  if (sequenceIn_>=0) {
5337    lowerIn_ = lower_[sequenceIn_];
5338    upperIn_ = upper_[sequenceIn_];
5339    valueIn_ = solution_[sequenceIn_];
5340    dualIn_ = dj_[sequenceIn_];
5341
5342    if (alpha_<0.0) {
5343      // as if from upper bound
5344      directionIn_=-1;
5345      upperIn_=valueIn_;
5346    } else {
5347      // as if from lower bound
5348      directionIn_=1;
5349      lowerIn_=valueIn_;
5350    }
5351  }
5352}
5353/*
5354   Row array has row part of pivot row
5355   Column array has column part.
5356   This is used in cleanup
5357*/
5358void
5359ClpSimplexDual::checkPossibleCleanup(CoinIndexedVector * rowArray,
5360                                        CoinIndexedVector * columnArray,
5361                                        double acceptablePivot)
5362{
5363  double * work;
5364  int number;
5365  int * which;
5366  int iSection;
5367
5368  double tolerance = dualTolerance_*1.001;
5369
5370  double thetaDown = 1.0e31;
5371  double thetaUp = 1.0e31;
5372  double bestAlphaDown = acceptablePivot*10.0;
5373  double bestAlphaUp = acceptablePivot*10.0;
5374  int sequenceDown =-1;
5375  int sequenceUp = -1;
5376
5377  double djSlack = dj_[pivotRow_];
5378  if (getRowStatus(pivotRow_)==basic)
5379    djSlack=COIN_DBL_MAX;
5380  if (fabs(djSlack)<tolerance)
5381    djSlack=0.0;
5382  int addSequence;
5383
5384  double alphaUp=0.0;
5385  double alphaDown=0.0;
5386  for (iSection=0;iSection<2;iSection++) {
5387
5388    int i;
5389    if (!iSection) {
5390      work = rowArray->denseVector();
5391      number = rowArray->getNumElements();
5392      which = rowArray->getIndices();
5393      addSequence = numberColumns_;
5394    } else {
5395      work = columnArray->denseVector();
5396      number = columnArray->getNumElements();
5397      which = columnArray->getIndices();
5398      addSequence = 0;
5399    }
5400   
5401    for (i=0;i<number;i++) {
5402      int iSequence = which[i];
5403      int iSequence2 = iSequence + addSequence;
5404      double alpha;
5405      double oldValue;
5406      double value;
5407
5408      switch(getStatus(iSequence2)) {
5409         
5410      case basic: 
5411        break;
5412      case ClpSimplex::isFixed:
5413        alpha = work[i];
5414        if (addSequence) {
5415          printf("possible - pivot row %d this %d\n",pivotRow_,iSequence);
5416          oldValue = dj_[iSequence2];
5417          if (alpha<=-acceptablePivot) {
5418            // might do other way
5419            value = oldValue+thetaUp*alpha;
5420            if (value<tolerance) {
5421              if (value<-tolerance||fabs(alpha)>bestAlphaUp) {
5422                thetaUp = -oldValue/alpha;
5423                bestAlphaUp = fabs(alpha);
5424                sequenceUp = iSequence2;
5425                alphaUp = alpha;
5426              }
5427            }
5428          } else if (alpha>=acceptablePivot) {
5429            // might do this way
5430            value = oldValue-thetaDown*alpha;
5431            if (value<tolerance) {
5432              if (value<-tolerance||fabs(alpha)>bestAlphaDown) {
5433                thetaDown = oldValue/alpha;
5434                bestAlphaDown = fabs(alpha);
5435                sequenceDown = iSequence2;
5436                alphaDown = alpha;
5437              }
5438            }
5439          }
5440        }
5441        break;
5442      case isFree:
5443      case superBasic:
5444        alpha = work[i];
5445        // dj must be effectively zero as dual feasible
5446        if (fabs(alpha)>bestAlphaUp) {
5447          thetaDown = 0.0;
5448          thetaUp = 0.0;
5449          bestAlphaDown = fabs(alpha);
5450          bestAlphaUp = bestAlphaUp;
5451          sequenceDown =iSequence2;
5452          sequenceUp = sequenceDown;
5453          alphaUp = alpha;
5454          alphaDown = alpha;
5455        }
5456        break;
5457      case atUpperBound:
5458        alpha = work[i];
5459        oldValue = dj_[iSequence2];
5460        if (alpha>=acceptablePivot) {
5461          // might do other way
5462          value = oldValue+thetaUp*alpha;
5463          if (value>-tolerance) {
5464            if (value>tolerance||fabs(alpha)>bestAlphaUp) {
5465              thetaUp = -oldValue/alpha;
5466              bestAlphaUp = fabs(alpha);
5467              sequenceUp = iSequence2;
5468              alphaUp = alpha;
5469            }
5470          }
5471        } else if (alpha<=-acceptablePivot) {
5472          // might do this way
5473          value = oldValue-thetaDown*alpha;
5474          if (value>-tolerance) {
5475            if (value>tolerance||fabs(alpha)>bestAlphaDown) {
5476              thetaDown = oldValue/alpha;
5477              bestAlphaDown = fabs(alpha);
5478              sequenceDown = iSequence2;
5479              alphaDown = alpha;
5480            }
5481          }
5482        }
5483        break;
5484      case atLowerBound:
5485        alpha = work[i];
5486        oldValue = dj_[iSequence2];
5487        if (alpha<=-acceptablePivot) {
5488          // might do other way
5489          value = oldValue+thetaUp*alpha;
5490          if (value<tolerance) {
5491            if (value<-tolerance||fabs(alpha)>bestAlphaUp) {
5492              thetaUp = -oldValue/alpha;
5493              bestAlphaUp = fabs(alpha);
5494              sequenceUp = iSequence2;
5495              alphaUp = alpha;
5496            }
5497          }
5498        } else if (alpha>=acceptablePivot) {
5499          // might do this way
5500          value = oldValue-thetaDown*alpha;
5501          if (value<tolerance) {
5502            if (value<-tolerance||fabs(alpha)>bestAlphaDown) {
5503              thetaDown = oldValue/alpha;
5504              bestAlphaDown = fabs(alpha);
5505              sequenceDown = iSequence2;
5506              alphaDown = alpha;
5507            }
5508          }
5509        }
5510        break;
5511      }
5512    }
5513  }
5514  thetaUp *= -1.0;
5515  // largest
5516  if (bestAlphaDown<bestAlphaUp) 
5517    sequenceDown =-1;
5518  else
5519    sequenceUp=-1;
5520
5521  sequenceIn_=-1;
5522 
5523  if (sequenceDown>=0) {
5524    theta_ = thetaDown;
5525    sequenceIn_ = sequenceDown;
5526    alpha_ = alphaDown;
5527#ifdef CLP_DEBUG
5528    if ((handler_->logLevel()&32))
5529      printf("predicted way - dirout %d, theta %g\n",
5530             directionOut_,theta_);
5531#endif
5532  } else if (sequenceUp>=0) {
5533    theta_ = thetaUp;
5534    sequenceIn_ = sequenceUp;
5535    alpha_ = alphaUp;
5536#ifdef CLP_DEBUG
5537    if ((handler_->logLevel()&32))
5538      printf("opposite way - dirout %d,theta %g\n",
5539             directionOut_,theta_);
5540#endif
5541  }
5542  if (sequenceIn_>=0) {
5543    lowerIn_ = lower_[sequenceIn_];
5544    upperIn_ = upper_[sequenceIn_];
5545    valueIn_ = solution_[sequenceIn_];
5546    dualIn_ = dj_[sequenceIn_];
5547
5548    if (alpha_<0.0) {
5549      // as if from upper bound
5550      directionIn_=-1;
5551      upperIn_=valueIn_;
5552    } else {
5553      // as if from lower bound
5554      directionIn_=1;
5555      lowerIn_=valueIn_;
5556    }
5557  }
5558}
5559/*
5560   This sees if we can move duals in dual values pass.
5561   This is done before any pivoting
5562*/
5563void ClpSimplexDual::doEasyOnesInValuesPass(double * dj)
5564{
5565  // Get column copy
5566  CoinPackedMatrix * columnCopy = matrix();
5567  // Get a row copy in standard format
5568  CoinPackedMatrix copy;
5569  copy.reverseOrderedCopyOf(*columnCopy);
5570  // get matrix data pointers
5571  const int * column = copy.getIndices();
5572  const CoinBigIndex * rowStart = copy.getVectorStarts();
5573  const int * rowLength = copy.getVectorLengths(); 
5574  const double * elementByRow = copy.getElements();
5575  double tolerance = dualTolerance_*1.001;
5576
5577  int iRow;
5578#ifdef CLP_DEBUG
5579  {
5580    double value5=0.0;
5581    int i;
5582    for (i=0;i<numberRows_+numberColumns_;i++) {
5583      if (dj[i]<-1.0e-6)
5584        value5 += dj[i]*upper_[i];
5585      else if (dj[i] >1.0e-6)
5586        value5 += dj[i]*lower_[i];
5587    }
5588    printf("Values objective Value before %g\n",value5);
5589  }
5590#endif
5591  // for scaled row
5592  double * scaled=NULL;
5593  if (rowScale_)
5594    scaled = new double[numberColumns_];
5595  for (iRow=0;iRow<numberRows_;iRow++) {
5596
5597    int iSequence = iRow + numberColumns_;
5598    double djBasic = dj[iSequence];
5599    if (getRowStatus(iRow)==basic&&fabs(djBasic)>tolerance) {
5600
5601      double changeUp ;
5602      // always -1 in pivot row
5603      if (djBasic>0.0) {
5604        // basic at lower bound
5605        changeUp = -lower_[iSequence];
5606      } else {
5607        // basic at upper bound
5608        changeUp = upper_[iSequence];
5609      }
5610      bool canMove=true;
5611      int i;
5612      const double * thisElements = elementByRow + rowStart[iRow]; 
5613      const int * thisIndices = column+rowStart[iRow];
5614      if (rowScale_) {
5615        assert (!auxiliaryModel_);
5616        // scale row
5617        double scale = rowScale_[iRow];
5618        for (i = 0;i<rowLength[iRow];i++) {
5619          int iColumn = thisIndices[i];
5620          double alpha = thisElements[i];
5621          scaled[i] = scale*alpha*columnScale_[iColumn];
5622        }
5623        thisElements = scaled;
5624      }
5625      for (i = 0;i<rowLength[iRow];i++) {
5626        int iColumn = thisIndices[i];
5627        double alpha = thisElements[i];
5628        double oldValue = dj[iColumn];;
5629        double value;
5630
5631        switch(getStatus(iColumn)) {
5632         
5633        case basic:
5634          if (dj[iColumn]<-tolerance&&
5635              fabs(solution_[iColumn]-upper_[iColumn])<1.0e-8) {
5636            // at ub
5637            changeUp += alpha*upper_[iColumn];
5638            // might do other way
5639            value = oldValue+djBasic*alpha;
5640            if (value>tolerance) 
5641              canMove=false;
5642          } else if (dj[iColumn]>tolerance&&
5643              fabs(solution_[iColumn]-lower_[iColumn])<1.0e-8) {
5644            changeUp += alpha*lower_[iColumn];
5645            // might do other way
5646            value = oldValue+djBasic*alpha;
5647            if (value<-tolerance)
5648              canMove=false;
5649          } else {
5650            canMove=false;
5651          }
5652          break;
5653        case ClpSimplex::isFixed:
5654          changeUp += alpha*upper_[iColumn];
5655          break;
5656        case isFree:
5657        case superBasic:
5658          canMove=false;
5659        break;
5660      case atUpperBound:
5661        changeUp += alpha*upper_[iColumn];
5662        // might do other way
5663        value = oldValue+djBasic*alpha;
5664        if (value>tolerance) 
5665          canMove=false;
5666        break;
5667      case atLowerBound:
5668        changeUp += alpha*lower_[iColumn];
5669        // might do other way
5670        value = oldValue+djBasic*alpha;
5671        if (value<-tolerance)
5672          canMove=false;
5673        break;
5674        }
5675      }
5676      if (canMove) {
5677        if (changeUp*djBasic>1.0e-12||fabs(changeUp)<1.0e-8) {
5678          // move
5679          for (i = 0;i<rowLength[iRow];i++) {
5680            int iColumn = thisIndices[i];
5681            double alpha = thisElements[i];
5682            dj[iColumn] += djBasic * alpha;
5683          }
5684          dj[iSequence]=0.0;
5685#ifdef CLP_DEBUG
5686          {
5687            double value5=0.0;
5688            int i;
5689            for (i=0;i<numberRows_+numberColumns_;i++) {
5690              if (dj[i]<-1.0e-6)
5691                value5 += dj[i]*upper_[i];
5692              else if (dj[i] >1.0e-6)
5693                value5 += dj[i]*lower_[i];
5694            }
5695            printf("Values objective Value after row %d old dj %g %g\n",
5696                   iRow,djBasic,value5);
5697          }
5698#endif
5699        }
5700      }
5701    }     
5702  }
5703  delete [] scaled;
5704}
5705int
5706ClpSimplexDual::nextSuperBasic()
5707{
5708  if (firstFree_>=0) {
5709    int returnValue=firstFree_;
5710    int iColumn=firstFree_+1;
5711    for (;iColumn<numberRows_+numberColumns_;iColumn++) {
5712      if (getStatus(iColumn)==isFree) 
5713        if (fabs(dj_[iColumn])>1.0e2*dualTolerance_) 
5714          break;
5715    }
5716    firstFree_ = iColumn;
5717    if (firstFree_==numberRows_+numberColumns_)
5718      firstFree_=-1;
5719    return returnValue;
5720  } else {
5721    return -1;
5722  }
5723}
5724void 
5725ClpSimplexDual::resetFakeBounds()
5726{
5727  double dummyChangeCost=0.0;
5728  changeBounds(true,rowArray_[2],dummyChangeCost);
5729  // throw away change
5730  for (int i=0;i<4;i++) 
5731    rowArray_[i]->clear();
5732}
Note: See TracBrowser for help on using the repository browser.