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

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

minor stuff

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 173.7 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        if (numberDualInfeasibilities_) {
3962          if (numberPrimalInfeasibilities_||numberPivots)
3963            problemStatus_=-1; // carry on as normal
3964          else
3965            problemStatus_=10; // try primal
3966        } else if (situationChanged==2) {
3967          problemStatus_=-1; // carry on as normal
3968        }
3969        situationChanged=0;
3970      } else {
3971        // iterate
3972        if (cleanDuals!=2) 
3973          problemStatus_=-1;
3974        else 
3975          problemStatus_ = 10; // try primal
3976      }
3977    }
3978  }
3979  if (type==0||type==1) {
3980    if (!type) {
3981      // create save arrays
3982      delete [] saveStatus_;
3983      delete [] savedSolution_;
3984      saveStatus_ = new unsigned char [numberRows_+numberColumns_];
3985      savedSolution_ = new double [numberRows_+numberColumns_];
3986    }
3987    // save arrays
3988    CoinMemcpyN(status_,numberColumns_+numberRows_,saveStatus_);
3989    CoinMemcpyN(rowActivityWork_,
3990           numberRows_,savedSolution_+numberColumns_);
3991    CoinMemcpyN(columnActivityWork_,numberColumns_,savedSolution_);
3992    // save extra stuff
3993    int dummy;
3994    matrix_->generalExpanded(this,5,dummy);
3995  }
3996
3997  // restore weights (if saved) - also recompute infeasibility list
3998  if (tentativeStatus>-3) 
3999    dualRowPivot_->saveWeights(this,(type <2) ? 2 : 4);
4000  else
4001    dualRowPivot_->saveWeights(this,3);
4002  // unflag all variables (we may want to wait a bit?)
4003  if ((tentativeStatus!=-2&&tentativeStatus!=-1)&&unflagVariables) {
4004    int iRow;
4005    int numberFlagged=0;
4006    for (iRow=0;iRow<numberRows_;iRow++) {
4007      int iPivot=pivotVariable_[iRow];
4008      if (flagged(iPivot)) {
4009        numberFlagged++;
4010        clearFlagged(iPivot);
4011      }
4012    }
4013#if 0
4014    if (numberFlagged) {
4015      printf("unflagging %d variables - status %d ninf %d nopt %d\n",numberFlagged,tentativeStatus,
4016             numberPrimalInfeasibilities_,
4017             numberTimesOptimal_);
4018    }
4019#endif
4020    unflagVariables = numberFlagged>0;
4021    if (numberFlagged&&!numberPivots) {
4022      /* looks like trouble as we have not done any iterations.
4023         Try changing pivot tolerance then give it a few goes and give up */
4024      if (factorization_->pivotTolerance()<0.9) {
4025        factorization_->pivotTolerance(0.99);
4026      } else if (numberTimesOptimal_<10) {
4027        numberTimesOptimal_++;
4028      } else {
4029        unflagVariables=false;
4030        changeMade_=false;
4031        secondaryStatus_ = 1; // and say probably infeasible
4032      }
4033    }
4034  }
4035  // see if cutoff reached
4036  double limit = 0.0;
4037  getDblParam(ClpDualObjectiveLimit, limit);
4038#if 0
4039  if(fabs(limit)<1.0e30&&objectiveValue()*optimizationDirection_>
4040     limit+100.0) {
4041    printf("lim %g obj %g %g - wo perturb %g sum dual %g\n",
4042           limit,objectiveValue_,objectiveValue(),computeInternalObjectiveValue(),sumDualInfeasibilities_);
4043  }
4044#endif
4045  if(fabs(limit)<1.0e30&&objectiveValue()*optimizationDirection_>
4046     limit&&!numberAtFakeBound()) {
4047    bool looksInfeasible = !numberDualInfeasibilities_;
4048    if (objectiveValue()*optimizationDirection_>limit+fabs(0.1*limit)+1.0e2*sumDualInfeasibilities_+1.0e4&&
4049        sumDualInfeasibilities_<largestDualError_&&numberIterations_>0.5*numberRows_+1000)
4050      looksInfeasible=true;
4051    if (looksInfeasible) {
4052      // Even if not perturbed internal costs may have changed
4053      if (true||perturbation_==101) {
4054        // be careful
4055        if (numberIterations_) {
4056          if(computeInternalObjectiveValue()>limit) {
4057            problemStatus_=1;
4058            secondaryStatus_ = 1; // and say was on cutoff
4059          }
4060        }
4061      } else {
4062        problemStatus_=1;
4063        secondaryStatus_ = 1; // and say was on cutoff
4064      }
4065    }
4066  }
4067  // If we are in trouble and in branch and bound give up
4068  if ((specialOptions_&1024)!=0) {
4069    int looksBad=0;
4070    if (largestPrimalError_*largestDualError_>1.0e2) {
4071      looksBad=1;
4072    } else if (largestPrimalError_>1.0e-2
4073        &&objectiveValue_>CoinMin(1.0e15,1.0e3*limit)) {
4074      looksBad=2;
4075    }
4076    if (looksBad) {
4077      if (factorization_->pivotTolerance()<0.9) {
4078        // up tolerance
4079        factorization_->pivotTolerance(CoinMin(factorization_->pivotTolerance()*1.05+0.02,0.91));
4080      } else if (numberIterations_>10000) {
4081        if (handler_->logLevel()>2)
4082          printf("bad dual - saying infeasible %d\n",looksBad);
4083        problemStatus_=1;
4084        secondaryStatus_ = 1; // and say was on cutoff
4085      } else if (largestPrimalError_>1.0e5) {
4086        {
4087          int iBigB=-1;
4088          double bigB=0.0;
4089          int iBigN=-1;
4090          double bigN=0.0;
4091          for (int i=0;i<numberRows_+numberColumns_;i++) {
4092            double value = fabs(solution_[i]);
4093            if (getStatus(i)==basic) {
4094              if (value>bigB) {
4095                bigB= value;
4096                iBigB=i;
4097              }
4098            } else {
4099              if (value>bigN) {
4100                bigN= value;
4101                iBigN=i;
4102              }
4103            }
4104          }
4105          if (bigB>1.0e8||bigN>1.0e8) {
4106            if (handler_->logLevel()>0)
4107              printf("it %d - basic %d %g, nonbasic %d %g\n",
4108                     numberIterations_,iBigB,bigB,iBigN,bigN);
4109          }
4110        }
4111        if (handler_->logLevel()>2)
4112          printf("bad dual - going to primal %d %g\n",looksBad,largestPrimalError_);
4113        allSlackBasis(true);
4114        problemStatus_=10;
4115      }
4116    }
4117  }
4118  if (problemStatus_<0&&!changeMade_) {
4119    problemStatus_=4; // unknown
4120  }
4121  lastGoodIteration_ = numberIterations_;
4122  if (problemStatus_<0) {
4123    sumDualInfeasibilities_=realDualInfeasibilities; // back to say be careful
4124    if (sumDualInfeasibilities_)
4125      numberDualInfeasibilities_=1;
4126  }
4127#if 1
4128  double thisObj = progress_->lastObjective(0);
4129  double lastObj = progress_->lastObjective(1);
4130  if (lastObj>thisObj+1.0e-4*CoinMax(fabs(thisObj),fabs(lastObj))+1.0e-4
4131      &&givenDuals==NULL&&firstFree_<0) {
4132    int maxFactor = factorization_->maximumPivots();
4133    if (maxFactor>10) {
4134      if (forceFactorization_<0)
4135        forceFactorization_= maxFactor;
4136      forceFactorization_ = CoinMax(1,(forceFactorization_>>1));
4137      //printf("Reducing factorization frequency\n");
4138    } 
4139  }
4140#endif
4141  // Allow matrices to be sorted etc
4142  int fake=-999; // signal sort
4143  matrix_->correctSequence(this,fake,fake);
4144}
4145/* While updateDualsInDual sees what effect is of flip
4146   this does actual flipping.
4147   If change >0.0 then value in array >0.0 => from lower to upper
4148*/
4149void 
4150ClpSimplexDual::flipBounds(CoinIndexedVector * rowArray,
4151                  CoinIndexedVector * columnArray,
4152                  double change)
4153{
4154  int number;
4155  int * which;
4156 
4157  int iSection;
4158
4159  for (iSection=0;iSection<2;iSection++) {
4160    int i;
4161    double * solution = solutionRegion(iSection);
4162    double * lower = lowerRegion(iSection);
4163    double * upper = upperRegion(iSection);
4164    int addSequence;
4165    if (!iSection) {
4166      number = rowArray->getNumElements();
4167      which = rowArray->getIndices();
4168      addSequence = numberColumns_;
4169    } else {
4170      number = columnArray->getNumElements();
4171      which = columnArray->getIndices();
4172      addSequence = 0;
4173    }
4174   
4175    for (i=0;i<number;i++) {
4176      int iSequence = which[i];
4177      Status status = getStatus(iSequence+addSequence);
4178
4179      switch(status) {
4180
4181      case basic:
4182      case isFree:
4183      case superBasic:
4184      case ClpSimplex::isFixed:
4185        break;
4186      case atUpperBound:
4187        // to lower bound
4188        setStatus(iSequence+addSequence,atLowerBound);
4189        solution[iSequence] = lower[iSequence];
4190        break;
4191      case atLowerBound:
4192        // to upper bound
4193        setStatus(iSequence+addSequence,atUpperBound);
4194        solution[iSequence] = upper[iSequence];
4195        break;
4196      }
4197    }
4198  }
4199  rowArray->setNumElements(0);
4200  columnArray->setNumElements(0);
4201}
4202// Restores bound to original bound
4203void 
4204ClpSimplexDual::originalBound( int iSequence)
4205{
4206  if (getFakeBound(iSequence)!=noFake) {
4207    numberFake_--;;
4208    setFakeBound(iSequence,noFake);
4209    if (auxiliaryModel_) {
4210      // just copy back
4211      lower_[iSequence]=auxiliaryModel_->lowerRegion()[iSequence+numberRows_+numberColumns_];
4212      upper_[iSequence]=auxiliaryModel_->upperRegion()[iSequence+numberRows_+numberColumns_];
4213      return;
4214    }
4215    if (iSequence>=numberColumns_) {
4216      // rows
4217      int iRow = iSequence-numberColumns_;
4218      rowLowerWork_[iRow]=rowLower_[iRow];
4219      rowUpperWork_[iRow]=rowUpper_[iRow];
4220      if (rowScale_) {
4221        if (rowLowerWork_[iRow]>-1.0e50)
4222          rowLowerWork_[iRow] *= rowScale_[iRow]*rhsScale_;
4223        if (rowUpperWork_[iRow]<1.0e50)
4224          rowUpperWork_[iRow] *= rowScale_[iRow]*rhsScale_;
4225      } else if (rhsScale_!=1.0) {
4226        if (rowLowerWork_[iRow]>-1.0e50)
4227          rowLowerWork_[iRow] *= rhsScale_;
4228        if (rowUpperWork_[iRow]<1.0e50)
4229          rowUpperWork_[iRow] *= rhsScale_;
4230      }
4231    } else {
4232      // columns
4233      columnLowerWork_[iSequence]=columnLower_[iSequence];
4234      columnUpperWork_[iSequence]=columnUpper_[iSequence];
4235      if (rowScale_) {
4236        double multiplier = 1.0/columnScale_[iSequence];
4237        if (columnLowerWork_[iSequence]>-1.0e50)
4238          columnLowerWork_[iSequence] *= multiplier*rhsScale_;
4239        if (columnUpperWork_[iSequence]<1.0e50)
4240          columnUpperWork_[iSequence] *= multiplier*rhsScale_;
4241      } else if (rhsScale_!=1.0) {
4242        if (columnLowerWork_[iSequence]>-1.0e50)
4243          columnLowerWork_[iSequence] *= rhsScale_;
4244        if (columnUpperWork_[iSequence]<1.0e50)
4245          columnUpperWork_[iSequence] *= rhsScale_;
4246      }
4247    }
4248  }
4249}
4250/* As changeBounds but just changes new bounds for a single variable.
4251   Returns true if change */
4252bool 
4253ClpSimplexDual::changeBound( int iSequence)
4254{
4255  // old values
4256  double oldLower=lower_[iSequence];
4257  double oldUpper=upper_[iSequence];
4258  double value=solution_[iSequence];
4259  bool modified=false;
4260  originalBound(iSequence);
4261  // original values
4262  double lowerValue=lower_[iSequence];
4263  double upperValue=upper_[iSequence];
4264  // back to altered values
4265  lower_[iSequence] = oldLower;
4266  upper_[iSequence] = oldUpper;
4267  assert (getFakeBound(iSequence)==noFake);
4268  //if (getFakeBound(iSequence)!=noFake)
4269  //numberFake_--;;
4270  if (value==oldLower) {
4271    if (upperValue > oldLower + dualBound_) {
4272      upper_[iSequence]=oldLower+dualBound_;
4273      setFakeBound(iSequence,upperFake);
4274      modified=true;
4275      numberFake_++;
4276    }
4277  } else if (value==oldUpper) {
4278    if (lowerValue < oldUpper - dualBound_) {
4279      lower_[iSequence]=oldUpper-dualBound_;
4280      setFakeBound(iSequence,lowerFake);
4281      modified=true;
4282      numberFake_++;
4283    }
4284  } else {
4285    assert(value==oldLower||value==oldUpper);
4286  }
4287  return modified;
4288}
4289// Perturbs problem
4290void 
4291ClpSimplexDual::perturb()
4292{
4293  if (perturbation_>100)
4294    return; //perturbed already
4295  if (perturbation_==100)
4296    perturbation_=50; // treat as normal
4297  int savePerturbation = perturbation_;
4298  bool modifyRowCosts=false;
4299  // dual perturbation
4300  double perturbation=1.0e-20;
4301  // maximum fraction of cost to perturb
4302  double maximumFraction = 1.0e-5;
4303  double constantPerturbation = 100.0*dualTolerance_;
4304  int maxLength=0;
4305  int minLength=numberRows_;
4306  double averageCost = 0.0;
4307  // look at element range
4308  double smallestNegative;
4309  double largestNegative;
4310  double smallestPositive;
4311  double largestPositive;
4312  matrix_->rangeOfElements(smallestNegative, largestNegative,
4313                           smallestPositive, largestPositive);
4314  smallestPositive = CoinMin(fabs(smallestNegative),smallestPositive);
4315  largestPositive = CoinMax(fabs(largestNegative),largestPositive);
4316  double elementRatio = largestPositive/smallestPositive;
4317  int numberNonZero=0;
4318  if (!numberIterations_&&perturbation_==50) {
4319    // See if we need to perturb
4320    double * sort = new double[numberColumns_];
4321    // Use objective BEFORE scaling
4322    const double * obj = objective();
4323    int i;
4324    for (i=0;i<numberColumns_;i++) {
4325      double value = fabs(obj[i]);
4326      sort[i]=value;
4327      averageCost += value;
4328      if (value)
4329        numberNonZero++;
4330    }
4331    if (numberNonZero)
4332      averageCost /= (double) numberNonZero;
4333    else
4334      averageCost = 1.0;
4335    std::sort(sort,sort+numberColumns_);
4336    int number=1;
4337    double last = sort[0];
4338    for (i=1;i<numberColumns_;i++) {
4339      if (last!=sort[i])
4340        number++;
4341      last=sort[i];
4342    }
4343    delete [] sort;
4344#if 0
4345    printf("nnz %d percent %d",number,(number*100)/numberColumns_);
4346    if (number*4>numberColumns_)
4347      printf(" - Would not perturb\n");
4348    else
4349      printf(" - Would perturb\n");
4350    //exit(0);
4351#endif
4352    //printf("ratio number diff costs %g, element ratio %g\n",((double)number)/((double) numberColumns_),
4353    //                                                                elementRatio);
4354    //number=0;
4355    if (number*4>numberColumns_||elementRatio>1.0e12) {
4356      perturbation_=100;
4357      return; // good enough
4358    }
4359  }
4360  int iColumn;
4361  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4362    if (columnLowerWork_[iColumn]<columnUpperWork_[iColumn]) {
4363      int length = matrix_->getVectorLength(iColumn);
4364      if (length>2) {
4365        maxLength = CoinMax(maxLength,length);
4366        minLength = CoinMin(minLength,length);
4367      }
4368    }
4369  }
4370  // If > 70 then do rows
4371  if (perturbation_>=70) {
4372    modifyRowCosts=true;
4373    perturbation_ -= 20;
4374    printf("Row costs modified, ");
4375  }
4376  bool uniformChange=false;
4377  if (perturbation_>50) {
4378    // Experiment
4379    // maximumFraction could be 1.0e-10 to 1.0
4380    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};
4381    maximumFraction = m[CoinMin(perturbation_-51,10)];
4382  }
4383  int iRow;
4384  double smallestNonZero=1.0e100;
4385  numberNonZero=0;
4386  if (perturbation_>=50) {
4387    perturbation = 1.0e-8;
4388    bool allSame=true;
4389    double lastValue=0.0;
4390    for (iRow=0;iRow<numberRows_;iRow++) {
4391      double lo = rowLowerWork_[iRow];
4392      double up = rowUpperWork_[iRow];
4393      if (lo<up) {
4394        double value = fabs(rowObjectiveWork_[iRow]);
4395        perturbation = CoinMax(perturbation,value);
4396        if (value) {
4397          modifyRowCosts=true;
4398          smallestNonZero = CoinMin(smallestNonZero,value);
4399        }
4400      } 
4401      if (lo&&lo>-1.0e10) {
4402        numberNonZero++;
4403        lo=fabs(lo);
4404        if (!lastValue) 
4405          lastValue=lo;
4406        else if (fabs(lo-lastValue)>1.0e-7)
4407          allSame=false;
4408      }
4409      if (up&&up<1.0e10) {
4410        numberNonZero++;
4411        up=fabs(up);
4412        if (!lastValue) 
4413          lastValue=up;
4414        else if (fabs(up-lastValue)>1.0e-7)
4415          allSame=false;
4416      }
4417    }
4418    double lastValue2=0.0;
4419    for (iColumn=0;iColumn<numberColumns_;iColumn++) { 
4420      double lo = columnLowerWork_[iColumn];
4421      double up = columnUpperWork_[iColumn];
4422      if (lo<up) {
4423        double value = 
4424          fabs(objectiveWork_[iColumn]);
4425        perturbation = CoinMax(perturbation,value);
4426        if (value) {
4427          smallestNonZero = CoinMin(smallestNonZero,value);
4428        }
4429      }
4430      if (lo&&lo>-1.0e10) {
4431        //numberNonZero++;
4432        lo=fabs(lo);
4433        if (!lastValue2) 
4434          lastValue2=lo;
4435        else if (fabs(lo-lastValue2)>1.0e-7)
4436          allSame=false;
4437      }
4438      if (up&&up<1.0e10) {
4439        //numberNonZero++;
4440        up=fabs(up);
4441        if (!lastValue2) 
4442          lastValue2=up;
4443        else if (fabs(up-lastValue2)>1.0e-7)
4444          allSame=false;
4445      }
4446    }
4447    if (allSame) {
4448      // Check elements
4449      double smallestNegative;
4450      double largestNegative;
4451      double smallestPositive;
4452      double largestPositive;
4453      matrix_->rangeOfElements(smallestNegative,largestNegative,
4454                               smallestPositive,largestPositive);
4455      if (smallestNegative==largestNegative&&
4456          smallestPositive==largestPositive) {
4457        // Really hit perturbation
4458        double adjust = CoinMin(100.0*maximumFraction,1.0e-3*CoinMax(lastValue,lastValue2));
4459        maximumFraction = CoinMax(adjust,maximumFraction);
4460      }
4461    }
4462    perturbation = CoinMin(perturbation,smallestNonZero/maximumFraction);
4463  } else {
4464    // user is in charge
4465    maximumFraction = 1.0e-1;
4466    // but some experiments
4467    if (perturbation_<=-900) {
4468      modifyRowCosts=true;
4469      perturbation_ += 1000;
4470      printf("Row costs modified, ");
4471    }
4472    if (perturbation_<=-10) {
4473      perturbation_ += 10; 
4474      maximumFraction = 1.0;
4475      if ((-perturbation_)%100>=10) {
4476        uniformChange=true;
4477        perturbation_+=20;
4478      }
4479      while (perturbation_<-10) {
4480        perturbation_ += 100;
4481        maximumFraction *= 1.0e-1;
4482      }
4483    }
4484    perturbation = pow(10.0,perturbation_);
4485  }
4486  double largestZero=0.0;
4487  double largest=0.0;
4488  double largestPerCent=0.0;
4489  // modify costs
4490  bool printOut=(handler_->logLevel()==63);
4491  printOut=false;
4492  if (modifyRowCosts) {
4493    for (iRow=0;iRow<numberRows_;iRow++) {
4494      if (rowLowerWork_[iRow]<rowUpperWork_[iRow]) {
4495        double value = perturbation;
4496        double currentValue = rowObjectiveWork_[iRow];
4497        value = CoinMin(value,maximumFraction*(fabs(currentValue)+1.0e-1*perturbation+1.0e-3));
4498        if (rowLowerWork_[iRow]>-largeValue_) {
4499          if (fabs(rowLowerWork_[iRow])<fabs(rowUpperWork_[iRow])) 
4500            value *= CoinDrand48();
4501          else
4502            value *= -CoinDrand48();
4503        } else if (rowUpperWork_[iRow]<largeValue_) {
4504          value *= -CoinDrand48();
4505        } else {
4506          value=0.0;
4507        }
4508        if (currentValue) {
4509          largest = CoinMax(largest,fabs(value));
4510          if (fabs(value)>fabs(currentValue)*largestPerCent)
4511            largestPerCent=fabs(value/currentValue);
4512        } else {
4513          largestZero=CoinMax(largestZero,fabs(value));
4514        }
4515        if (printOut)
4516          printf("row %d cost %g change %g\n",iRow,rowObjectiveWork_[iRow],value);
4517        rowObjectiveWork_[iRow] += value;
4518      }
4519    }
4520  }
4521  // 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};
4522  // 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};
4523  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};
4524  //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};
4525  double extraWeight=10.0;
4526  // Scale back if wanted
4527  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};
4528  if (constantPerturbation<99.0*dualTolerance_) {
4529    perturbation *= 0.1;
4530    extraWeight=0.5;
4531    memcpy(weight,weight2,sizeof(weight2));
4532  }
4533  // adjust weights if all columns long
4534  double factor=1.0;
4535  if (maxLength) {
4536    factor = 3.0/(double) minLength;
4537  }
4538  // Make variables with more elements more expensive
4539  const double m1 = 0.5;
4540  double smallestAllowed = CoinMin(1.0e-2*dualTolerance_,maximumFraction);
4541  //double largestAllowed = CoinMax(1.0e3*dualTolerance_,maximumFraction*10.0*averageCost);
4542  double largestAllowed = CoinMax(1.0e3*dualTolerance_,maximumFraction*averageCost);
4543  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4544    if (columnLowerWork_[iColumn]<columnUpperWork_[iColumn]&&getStatus(iColumn)!=basic) {
4545      double value = perturbation;
4546      double currentValue = objectiveWork_[iColumn];
4547      value = CoinMin(value,constantPerturbation+maximumFraction*(fabs(currentValue)+1.0e-1*perturbation+1.0e-8));
4548      //value = CoinMin(value,constantPerturbation;+maximumFraction*fabs(currentValue));
4549      double value2 = constantPerturbation+1.0e-1*smallestNonZero;
4550      if (uniformChange) {
4551        value = maximumFraction;
4552        value2=maximumFraction;
4553      }
4554      if (columnLowerWork_[iColumn]>-largeValue_) {
4555        if (fabs(columnLowerWork_[iColumn])<
4556            fabs(columnUpperWork_[iColumn])) {
4557          value *= (1.0-m1 +m1*CoinDrand48());
4558          value2 *= (1.0-m1 +m1*CoinDrand48());
4559        } else {
4560          value *= -(1.0-m1+m1*CoinDrand48());
4561          value2 *= -(1.0-m1+m1*CoinDrand48());
4562        }
4563      } else if (columnUpperWork_[iColumn]<largeValue_) {
4564        value *= -(1.0-m1+m1*CoinDrand48());
4565        value2 *= -(1.0-m1+m1*CoinDrand48());
4566      } else {
4567        value=0.0;
4568      }
4569      if (value) {
4570        int length = matrix_->getVectorLength(iColumn);
4571        if (length>3) {
4572          length = (int) ((double) length * factor);
4573          length = CoinMax(3,length);
4574        }
4575        double multiplier;
4576#if 1
4577        if (length<10)
4578          multiplier=weight[length];
4579        else
4580          multiplier = weight[10];
4581#else
4582        if (length<10)
4583          multiplier=weight[length];
4584        else
4585          multiplier = weight[10]+extraWeight*(length-10);
4586        multiplier *= 0.5;
4587#endif
4588        value *= multiplier;
4589        value = CoinMin(value,value2);
4590        if (savePerturbation<50||savePerturbation>60) {
4591          if (fabs(value)<=dualTolerance_)
4592            value=0.0;
4593        } else if (value) {
4594          // get in range
4595          if (fabs(value)<=smallestAllowed) {
4596            value *= 10.0;
4597            while (fabs(value)<=smallestAllowed) 
4598              value *= 10.0;
4599          } else if (fabs(value)>largestAllowed) {
4600            value *= 0.1;
4601            while (fabs(value)>largestAllowed) 
4602              value *= 0.1;
4603          }
4604        }
4605        if (currentValue) {
4606          largest = CoinMax(largest,fabs(value));
4607          if (fabs(value)>fabs(currentValue)*largestPerCent)
4608            largestPerCent=fabs(value/currentValue);
4609        } else {
4610          largestZero=CoinMax(largestZero,fabs(value));
4611        }
4612        // but negative if at ub
4613        if (getStatus(iColumn)==atUpperBound)
4614          value = -value;
4615        if (printOut)
4616          printf("col %d cost %g change %g\n",iColumn,objectiveWork_[iColumn],value);
4617        objectiveWork_[iColumn] += value;
4618      }
4619    }
4620  }
4621  handler_->message(CLP_SIMPLEX_PERTURB,messages_)
4622    <<100.0*maximumFraction<<perturbation<<largest<<100.0*largestPerCent<<largestZero
4623    <<CoinMessageEol;
4624  // and zero changes
4625  //int nTotal = numberRows_+numberColumns_;
4626  //CoinZeroN(cost_+nTotal,nTotal);
4627  // say perturbed
4628  perturbation_=101;
4629
4630}
4631/* For strong branching.  On input lower and upper are new bounds
4632   while on output they are change in objective function values
4633   (>1.0e50 infeasible).
4634   Return code is 0 if nothing interesting, -1 if infeasible both
4635   ways and +1 if infeasible one way (check values to see which one(s))
4636   Returns -2 if bad factorization
4637*/
4638int ClpSimplexDual::strongBranching(int numberVariables,const int * variables,
4639                                    double * newLower, double * newUpper,
4640                                    double ** outputSolution,
4641                                    int * outputStatus, int * outputIterations,
4642                                    bool stopOnFirstInfeasible,
4643                                    bool alwaysFinish,
4644                                    int startFinishOptions)
4645{
4646  int i;
4647  int returnCode=0;
4648  double saveObjectiveValue = objectiveValue_;
4649  algorithm_ = -1;
4650
4651  //scaling(false);
4652
4653  // put in standard form (and make row copy)
4654  // create modifiable copies of model rim and do optional scaling
4655  createRim(7+8+16+32,true,startFinishOptions);
4656
4657  // change newLower and newUpper if scaled
4658
4659  // Do initial factorization
4660  // and set certain stuff
4661  // We can either set increasing rows so ...IsBasic gives pivot row
4662  // or we can just increment iBasic one by one
4663  // for now let ...iBasic give pivot row
4664  int useFactorization=false;
4665  if ((startFinishOptions&2)!=0&&(whatsChanged_&(2+512))==2+512)
4666    useFactorization=true; // Keep factorization if possible
4667  // switch off factorization if bad
4668  if (pivotVariable_[0]<0)
4669    useFactorization=false;
4670  if (!useFactorization||factorization_->numberRows()!=numberRows_) {
4671    useFactorization = false;
4672    factorization_->increasingRows(2);
4673    // row activities have negative sign
4674    factorization_->slackValue(-1.0);
4675    factorization_->zeroTolerance(1.0e-13);
4676
4677    int factorizationStatus = internalFactorize(0);
4678    if (factorizationStatus<0) {
4679      // some error
4680      // we should either debug or ignore
4681#ifndef NDEBUG
4682      printf("***** ClpDual strong branching factorization error - debug\n");
4683#endif
4684      return -2;
4685    } else if (factorizationStatus&&factorizationStatus<=numberRows_) {
4686      handler_->message(CLP_SINGULARITIES,messages_)
4687        <<factorizationStatus
4688        <<CoinMessageEol;
4689    }
4690  }
4691  // save stuff
4692  ClpFactorization saveFactorization(*factorization_);
4693  // save basis and solution
4694  double * saveSolution = new double[numberRows_+numberColumns_];
4695  CoinMemcpyN(solution_,
4696         numberRows_+numberColumns_,saveSolution);
4697  unsigned char * saveStatus = 
4698    new unsigned char [numberRows_+numberColumns_];
4699  CoinMemcpyN(status_,numberColumns_+numberRows_,saveStatus);
4700  // save bounds as createRim makes clean copies
4701  double * saveLower = new double[numberRows_+numberColumns_];
4702  CoinMemcpyN(lower_,
4703         numberRows_+numberColumns_,saveLower);
4704  double * saveUpper = new double[numberRows_+numberColumns_];
4705  CoinMemcpyN(upper_,
4706         numberRows_+numberColumns_,saveUpper);
4707  double * saveObjective = new double[numberRows_+numberColumns_];
4708  CoinMemcpyN(cost_,
4709         numberRows_+numberColumns_,saveObjective);
4710  int * savePivot = new int [numberRows_];
4711  CoinMemcpyN(pivotVariable_, numberRows_,savePivot);
4712  // need to save/restore weights.
4713
4714  int iSolution = 0;
4715  for (i=0;i<numberVariables;i++) {
4716    int iColumn = variables[i];
4717    double objectiveChange;
4718    double saveBound;
4719   
4720    // try down
4721   
4722    saveBound = columnUpper_[iColumn];
4723    // external view - in case really getting optimal
4724    columnUpper_[iColumn] = newUpper[i];
4725    if (scalingFlag_<=0) 
4726      upper_[iColumn] = newUpper[i]*rhsScale_;
4727    else 
4728      upper_[iColumn] = (newUpper[i]/columnScale_[iColumn])*rhsScale_; // scale
4729    // Start of fast iterations
4730    int status = fastDual(alwaysFinish);
4731    CoinAssert (problemStatus_||objectiveValue_<1.0e50);
4732    // make sure plausible
4733    double obj = CoinMax(objectiveValue_,saveObjectiveValue);
4734    if (status&&problemStatus_!=3) {
4735      // not finished - might be optimal
4736      checkPrimalSolution(rowActivityWork_,columnActivityWork_);
4737      double limit = 0.0;
4738      getDblParam(ClpDualObjectiveLimit, limit);
4739      if (!numberPrimalInfeasibilities_&&obj<limit) { 
4740        problemStatus_=0;
4741      } 
4742      status=problemStatus_;
4743    }
4744    if (status||(problemStatus_==0&&!isDualObjectiveLimitReached())) {
4745      objectiveChange = obj-saveObjectiveValue;
4746    } else {
4747      objectiveChange = 1.0e100;
4748      status=1;
4749    }
4750    if (problemStatus_==3)
4751      status=2;
4752
4753    if (scalingFlag_<=0) {
4754      CoinMemcpyN(solution_,numberColumns_,outputSolution[iSolution]);
4755    } else {
4756      int j;
4757      double * sol = outputSolution[iSolution];
4758      for (j=0;j<numberColumns_;j++) 
4759        sol[j] = solution_[j]*columnScale_[j];
4760    }
4761    outputStatus[iSolution]=status;
4762    outputIterations[iSolution]=numberIterations_;
4763    iSolution++;
4764    // restore
4765    CoinMemcpyN(saveSolution,
4766            numberRows_+numberColumns_,solution_);
4767    CoinMemcpyN(saveStatus,numberColumns_+numberRows_,status_);
4768    CoinMemcpyN(saveLower,
4769            numberRows_+numberColumns_,lower_);
4770    CoinMemcpyN(saveUpper,
4771            numberRows_+numberColumns_,upper_);
4772    CoinMemcpyN(saveObjective,
4773            numberRows_+numberColumns_,cost_);
4774    columnUpper_[iColumn] = saveBound;
4775    CoinMemcpyN(savePivot, numberRows_,pivotVariable_);
4776    delete factorization_;
4777    factorization_ = new ClpFactorization(saveFactorization);
4778
4779    newUpper[i]=objectiveChange;
4780#ifdef CLP_DEBUG
4781    printf("down on %d costs %g\n",iColumn,objectiveChange);
4782#endif
4783         
4784    // try up
4785   
4786    saveBound = columnLower_[iColumn];
4787    // external view - in case really getting optimal
4788    columnLower_[iColumn] = newLower[i];
4789    if (scalingFlag_<=0) 
4790      lower_[iColumn] = newLower[i]*rhsScale_;
4791    else 
4792      lower_[iColumn] = (newLower[i]/columnScale_[iColumn])*rhsScale_; // scale
4793    // Start of fast iterations
4794    status = fastDual(alwaysFinish);
4795    // make sure plausible
4796    obj = CoinMax(objectiveValue_,saveObjectiveValue);
4797    if (status&&problemStatus_!=3) {
4798      // not finished - might be optimal
4799      checkPrimalSolution(rowActivityWork_,columnActivityWork_);
4800      double limit = 0.0;
4801      getDblParam(ClpDualObjectiveLimit, limit);
4802      if (!numberPrimalInfeasibilities_&&obj< limit) { 
4803        problemStatus_=0;
4804      } 
4805      status=problemStatus_;
4806    }
4807    if (status||(problemStatus_==0&&!isDualObjectiveLimitReached())) {
4808      objectiveChange = obj-saveObjectiveValue;
4809    } else {
4810      objectiveChange = 1.0e100;
4811      status=1;
4812    }
4813    if (problemStatus_==3)
4814      status=2;
4815    if (scalingFlag_<=0) {
4816      CoinMemcpyN(solution_,numberColumns_,outputSolution[iSolution]);
4817    } else {
4818      int j;
4819      double * sol = outputSolution[iSolution];
4820      for (j=0;j<numberColumns_;j++) 
4821        sol[j] = solution_[j]*columnScale_[j];
4822    }
4823    outputStatus[iSolution]=status;
4824    outputIterations[iSolution]=numberIterations_;
4825    iSolution++;
4826
4827    // restore
4828    CoinMemcpyN(saveSolution,
4829            numberRows_+numberColumns_,solution_);
4830    CoinMemcpyN(saveStatus,numberColumns_+numberRows_,status_);
4831    CoinMemcpyN(saveLower,
4832            numberRows_+numberColumns_,lower_);
4833    CoinMemcpyN(saveUpper,
4834            numberRows_+numberColumns_,upper_);
4835    CoinMemcpyN(saveObjective,
4836            numberRows_+numberColumns_,cost_);
4837    columnLower_[iColumn] = saveBound;
4838    CoinMemcpyN(savePivot, numberRows_,pivotVariable_);
4839    delete factorization_;
4840    factorization_ = new ClpFactorization(saveFactorization);
4841
4842    newLower[i]=objectiveChange;
4843#ifdef CLP_DEBUG
4844    printf("up on %d costs %g\n",iColumn,objectiveChange);
4845#endif
4846         
4847    /* Possibilities are:
4848       Both sides feasible - store
4849       Neither side feasible - set objective high and exit if desired
4850       One side feasible - change bounds and resolve
4851    */
4852    if (newUpper[i]<1.0e100) {
4853      if(newLower[i]<1.0e100) {
4854        // feasible - no action
4855      } else {
4856        // up feasible, down infeasible
4857        returnCode=1;
4858        if (stopOnFirstInfeasible)
4859          break;
4860      }
4861    } else {
4862      if(newLower[i]<1.0e100) {
4863        // down feasible, up infeasible
4864        returnCode=1;
4865        if (stopOnFirstInfeasible)
4866          break;
4867      } else {
4868        // neither side feasible
4869        returnCode=-1;
4870        break;
4871      }
4872    }
4873  }
4874  delete [] saveSolution;
4875  delete [] saveLower;
4876  delete [] saveUpper;
4877  delete [] saveObjective;
4878  delete [] saveStatus;
4879  delete [] savePivot;
4880  if ((startFinishOptions&1)==0) {
4881    deleteRim(1);
4882    whatsChanged_=0;
4883  } else {
4884    // Original factorization will have been put back by last loop
4885    //delete factorization_;
4886    //factorization_ = new ClpFactorization(saveFactorization);
4887    deleteRim(0);
4888    // mark all as current
4889    whatsChanged_ = 0xffff;
4890  }
4891  objectiveValue_ = saveObjectiveValue;
4892  return returnCode;
4893}
4894// treat no pivot as finished (unless interesting)
4895int ClpSimplexDual::fastDual(bool alwaysFinish)
4896{
4897  algorithm_ = -1;
4898  secondaryStatus_=0;
4899  // Say in fast dual
4900  specialOptions_ |= 16384;
4901  //handler_->setLogLevel(63);
4902  // save data
4903  ClpDataSave data = saveData();
4904  dualTolerance_=dblParam_[ClpDualTolerance];
4905  primalTolerance_=dblParam_[ClpPrimalTolerance];
4906
4907  // save dual bound
4908  double saveDualBound = dualBound_;
4909
4910  double objectiveChange;
4911  // for dual we will change bounds using dualBound_
4912  // for this we need clean basis so it is after factorize
4913  gutsOfSolution(NULL,NULL);
4914  numberFake_ =0; // Number of variables at fake bounds
4915  numberChanged_ =0; // Number of variables with changed costs
4916  changeBounds(true,NULL,objectiveChange);
4917
4918  problemStatus_ = -1;
4919  numberIterations_=0;
4920  factorization_->sparseThreshold(0);
4921  factorization_->goSparse();
4922
4923  int lastCleaned=0; // last time objective or bounds cleaned up
4924
4925  // number of times we have declared optimality
4926  numberTimesOptimal_=0;
4927
4928  // This says whether to restore things etc
4929  int factorType=0;
4930  /*
4931    Status of problem:
4932    0 - optimal
4933    1 - infeasible
4934    2 - unbounded
4935    -1 - iterating
4936    -2 - factorization wanted
4937    -3 - redo checking without factorization
4938    -4 - looks infeasible
4939
4940    BUT also from whileIterating return code is:
4941
4942   -1 iterations etc
4943   -2 inaccuracy
4944   -3 slight inaccuracy (and done iterations)
4945   +0 looks optimal (might be unbounded - but we will investigate)
4946   +1 looks infeasible
4947   +3 max iterations
4948
4949  */
4950
4951  int returnCode = 0;
4952
4953  int iRow,iColumn;
4954  while (problemStatus_<0) {
4955    // clear
4956    for (iRow=0;iRow<4;iRow++) {
4957      rowArray_[iRow]->clear();
4958    }   
4959   
4960    for (iColumn=0;iColumn<2;iColumn++) {
4961      columnArray_[iColumn]->clear();
4962    }   
4963
4964    // give matrix (and model costs and bounds a chance to be
4965    // refreshed (normally null)
4966    matrix_->refresh(this);
4967    // may factorize, checks if problem finished
4968    // should be able to speed this up on first time
4969    statusOfProblemInDual(lastCleaned,factorType,NULL,data,0);
4970
4971    // Say good factorization
4972    factorType=1;
4973
4974    // Do iterations
4975    if (problemStatus_<0) {
4976      double * givenPi=NULL;
4977      returnCode = whileIterating(givenPi,0);
4978      if ((!alwaysFinish&&returnCode<0)||returnCode==3) {
4979        if (returnCode!=3)
4980          assert (problemStatus_<0);
4981        returnCode=1;
4982        problemStatus_ = 3; 
4983        // can't say anything interesting - might as well return
4984#ifdef CLP_DEBUG
4985        printf("returning from fastDual after %d iterations with code %d\n",
4986               numberIterations_,returnCode);
4987#endif
4988        break;
4989      }
4990      returnCode=0;
4991    }
4992  }
4993
4994  // clear
4995  for (iRow=0;iRow<4;iRow++) {
4996    rowArray_[iRow]->clear();
4997  }   
4998 
4999  for (iColumn=0;iColumn<2;iColumn++) {
5000    columnArray_[iColumn]->clear();
5001  }   
5002  // Say not in fast dual
5003  specialOptions_ &= ~16384;
5004  assert(!numberFake_||((specialOptions_&(2048|4096))!=0&&dualBound_>1.0e8)
5005         ||returnCode||problemStatus_); // all bounds should be okay
5006  // Restore any saved stuff
5007  restoreData(data);
5008  dualBound_ = saveDualBound;
5009  return returnCode;
5010}
5011// This does first part of StrongBranching
5012ClpFactorization * 
5013ClpSimplexDual::setupForStrongBranching(char * arrays, int numberRows, int numberColumns)
5014{
5015  algorithm_ = -1;
5016  // put in standard form (and make row copy)
5017  // create modifiable copies of model rim and do optional scaling
5018  int startFinishOptions;
5019  /*  COIN_CLP_VETTED
5020      Looks safe for Cbc
5021  */
5022  bool safeOption = ((specialOptions_&COIN_CBC_USING_CLP)!=0);
5023  if((specialOptions_&4096)==0||(auxiliaryModel_&&!safeOption)) {
5024    startFinishOptions=0;
5025  } else {
5026    startFinishOptions=1+2+4;
5027  }
5028  createRim(7+8+16+32,true,startFinishOptions);
5029  // Do initial factorization
5030  // and set certain stuff
5031  // We can either set increasing rows so ...IsBasic gives pivot row
5032  // or we can just increment iBasic one by one
5033  // for now let ...iBasic give pivot row
5034  bool useFactorization=false;
5035  if ((startFinishOptions&2)!=0&&(whatsChanged_&(2+512))==2+512) {
5036    useFactorization=true; // Keep factorization if possible
5037    // switch off factorization if bad
5038    if (pivotVariable_[0]<0||factorization_->numberRows()!=numberRows_) 
5039      useFactorization=false;
5040  }
5041  if (!useFactorization) {
5042    factorization_->increasingRows(2);
5043    // row activities have negative sign
5044    factorization_->slackValue(-1.0);
5045    factorization_->zeroTolerance(1.0e-13);
5046
5047    int factorizationStatus = internalFactorize(0);
5048    if (factorizationStatus<0) {
5049      // some error
5050      // we should either debug or ignore
5051#ifndef NDEBUG
5052      printf("***** ClpDual strong branching factorization error - debug\n");
5053#endif
5054    } else if (factorizationStatus&&factorizationStatus<=numberRows_) {
5055      handler_->message(CLP_SINGULARITIES,messages_)
5056        <<factorizationStatus
5057        <<CoinMessageEol;
5058    }
5059  }
5060  double * arrayD = (double *) arrays;
5061  arrayD[0]=objectiveValue()*optimizationDirection_;
5062  double * saveSolution = arrayD+1;
5063  double * saveLower = saveSolution + (numberRows+numberColumns);
5064  double * saveUpper = saveLower + (numberRows+numberColumns);
5065  double * saveObjective = saveUpper + (numberRows+numberColumns);
5066  double * saveLowerOriginal = saveObjective + (numberRows+numberColumns);
5067  double * saveUpperOriginal = saveLowerOriginal + numberColumns;
5068  arrayD = saveUpperOriginal + numberColumns;
5069  int * savePivot = (int *) arrayD;
5070  int * whichRow = savePivot+numberRows;
5071  int * whichColumn = whichRow + 3*numberRows;
5072  int * arrayI = whichColumn + 2*numberColumns;
5073  unsigned char * saveStatus = (unsigned char *) (arrayI+1);
5074  // save stuff
5075  // save basis and solution
5076  CoinMemcpyN(solution_,
5077          numberRows_+numberColumns_,saveSolution);
5078  CoinMemcpyN(status_,numberColumns_+numberRows_,saveStatus);
5079  CoinMemcpyN(lower_,
5080          numberRows_+numberColumns_,saveLower);
5081  CoinMemcpyN(upper_,
5082          numberRows_+numberColumns_,saveUpper);
5083  CoinMemcpyN(cost_,
5084          numberRows_+numberColumns_,saveObjective);
5085  CoinMemcpyN(pivotVariable_, numberRows_,savePivot);
5086  return new ClpFactorization(*factorization_);
5087}
5088// This cleans up after strong branching
5089void 
5090ClpSimplexDual::cleanupAfterStrongBranching()
5091{
5092  int startFinishOptions;
5093  /*  COIN_CLP_VETTED
5094      Looks safe for Cbc
5095  */
5096  bool safeOption = ((specialOptions_&COIN_CBC_USING_CLP)!=0);
5097  if((specialOptions_&4096)==0||(auxiliaryModel_&&!safeOption)) {
5098    startFinishOptions=0;
5099  } else {
5100    startFinishOptions=1+2+4;
5101  }
5102  if ((startFinishOptions&1)==0) {
5103    deleteRim(1);
5104    whatsChanged_=0;
5105  } else {
5106    // Original factorization will have been put back by last loop
5107    //delete factorization_;
5108    //factorization_ = new ClpFactorization(saveFactorization);
5109    deleteRim(0);
5110    // mark all as current
5111    whatsChanged_ = 0xffff;
5112  }
5113}
5114/* Checks number of variables at fake bounds.  This is used by fastDual
5115   so can exit gracefully before end */
5116int 
5117ClpSimplexDual::numberAtFakeBound()
5118{
5119  int iSequence;
5120  int numberFake=0;
5121 
5122  for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
5123    FakeBound bound = getFakeBound(iSequence);
5124    switch(getStatus(iSequence)) {
5125
5126    case basic:
5127      break;
5128    case isFree:
5129    case superBasic:
5130    case ClpSimplex::isFixed:
5131      //setFakeBound (iSequence, noFake);
5132      break;
5133    case atUpperBound:
5134      if (bound==upperFake||bound==bothFake)
5135        numberFake++;
5136      break;
5137    case atLowerBound:
5138      if (bound==lowerFake||bound==bothFake)
5139        numberFake++;
5140      break;
5141    }
5142  }
5143  numberFake_ = numberFake;
5144  return numberFake;
5145}
5146/* Pivot out a variable and choose an incoing one.  Assumes dual
5147   feasible - will not go through a reduced cost. 
5148   Returns step length in theta
5149   Returns ray in ray_ (or NULL if no pivot)
5150   Return codes as before but -1 means no acceptable pivot
5151*/
5152int 
5153ClpSimplexDual::pivotResult()
5154{
5155  abort();
5156  return 0;
5157}
5158/*
5159   Row array has row part of pivot row
5160   Column array has column part.
5161   This is used in dual values pass
5162*/
5163void
5164ClpSimplexDual::checkPossibleValuesMove(CoinIndexedVector * rowArray,
5165                                        CoinIndexedVector * columnArray,
5166                                        double acceptablePivot)
5167{
5168  double * work;
5169  int number;
5170  int * which;
5171  int iSection;
5172
5173  double tolerance = dualTolerance_*1.001;
5174
5175  double thetaDown = 1.0e31;
5176  double changeDown ;
5177  double thetaUp = 1.0e31;
5178  double bestAlphaDown = acceptablePivot*0.99999;
5179  double bestAlphaUp = acceptablePivot*0.99999;
5180  int sequenceDown =-1;
5181  int sequenceUp = sequenceOut_;
5182
5183  double djBasic = dj_[sequenceOut_];
5184  if (djBasic>0.0) {
5185    // basic at lower bound so directionOut_ 1 and -1 in pivot row
5186    // dj will go to zero on other way
5187    thetaUp = djBasic;
5188    changeDown = -lower_[sequenceOut_];
5189  } else {
5190    // basic at upper bound so directionOut_ -1 and 1 in pivot row
5191    // dj will go to zero on other way
5192    thetaUp = -djBasic;
5193    changeDown = upper_[sequenceOut_];
5194  }
5195  bestAlphaUp = 1.0;
5196  int addSequence;
5197
5198  double alphaUp=0.0;
5199  double alphaDown=0.0;
5200
5201  for (iSection=0;iSection<2;iSection++) {
5202
5203    int i;
5204    if (!iSection) {
5205      work = rowArray->denseVector();
5206      number = rowArray->getNumElements();
5207      which = rowArray->getIndices();
5208      addSequence = numberColumns_;
5209    } else {
5210      work = columnArray->denseVector();
5211      number = columnArray->getNumElements();
5212      which = columnArray->getIndices();
5213      addSequence = 0;
5214    }
5215   
5216    for (i=0;i<number;i++) {
5217      int iSequence = which[i];
5218      int iSequence2 = iSequence + addSequence;
5219      double alpha;
5220      double oldValue;
5221      double value;
5222
5223      switch(getStatus(iSequence2)) {
5224         
5225      case basic: 
5226        break;
5227      case ClpSimplex::isFixed:
5228        alpha = work[i];
5229        changeDown += alpha*upper_[iSequence2];
5230        break;
5231      case isFree:
5232      case superBasic:
5233        alpha = work[i];
5234        // dj must be effectively zero as dual feasible
5235        if (fabs(alpha)>bestAlphaUp) {
5236          thetaDown = 0.0;
5237          thetaUp = 0.0;
5238          bestAlphaDown = fabs(alpha);
5239          bestAlphaUp = bestAlphaDown;
5240          sequenceDown =iSequence2;
5241          sequenceUp = sequenceDown;
5242          alphaUp = alpha;
5243          alphaDown = alpha;
5244        }
5245        break;
5246      case atUpperBound:
5247        alpha = work[i];
5248        oldValue = dj_[iSequence2];
5249        changeDown += alpha*upper_[iSequence2];
5250        if (alpha>=acceptablePivot) {
5251          // might do other way
5252          value = oldValue+thetaUp*alpha;
5253          if (value>-tolerance) {
5254            if (value>tolerance||fabs(alpha)>bestAlphaUp) {
5255              thetaUp = -oldValue/alpha;
5256              bestAlphaUp = fabs(alpha);
5257              sequenceUp = iSequence2;
5258              alphaUp = alpha;
5259            }
5260          }
5261        } else if (alpha<=-acceptablePivot) {
5262          // might do this way
5263          value = oldValue-thetaDown*alpha;
5264          if (value>-tolerance) {
5265            if (value>tolerance||fabs(alpha)>bestAlphaDown) {
5266              thetaDown = oldValue/alpha;
5267              bestAlphaDown = fabs(alpha);
5268              sequenceDown = iSequence2;
5269              alphaDown = alpha;
5270            }
5271          }
5272        }
5273        break;
5274      case atLowerBound:
5275        alpha = work[i];
5276        oldValue = dj_[iSequence2];
5277        changeDown += alpha*lower_[iSequence2];
5278        if (alpha<=-acceptablePivot) {
5279          // might do other way
5280          value = oldValue+thetaUp*alpha;
5281          if (value<tolerance) {
5282            if (value<-tolerance||fabs(alpha)>bestAlphaUp) {
5283              thetaUp = -oldValue/alpha;
5284              bestAlphaUp = fabs(alpha);
5285              sequenceUp = iSequence2;
5286              alphaUp = alpha;
5287            }
5288          }
5289        } else if (alpha>=acceptablePivot) {
5290          // might do this way
5291          value = oldValue-thetaDown*alpha;
5292          if (value<tolerance) {
5293            if (value<-tolerance||fabs(alpha)>bestAlphaDown) {
5294              thetaDown = oldValue/alpha;
5295              bestAlphaDown = fabs(alpha);
5296              sequenceDown = iSequence2;
5297              alphaDown = alpha;
5298            }
5299          }
5300        }
5301        break;
5302      }
5303    }
5304  }
5305  thetaUp *= -1.0;
5306  double changeUp = -thetaUp * changeDown;
5307  changeDown = -thetaDown*changeDown;
5308  if (CoinMax(fabs(thetaDown),fabs(thetaUp))<1.0e-8) {
5309    // largest
5310    if (fabs(alphaDown)<fabs(alphaUp)) {
5311      sequenceDown =-1;
5312    }
5313  }
5314  // choose
5315  sequenceIn_=-1;
5316  if (changeDown>changeUp&&sequenceDown>=0) {
5317    theta_ = thetaDown;
5318    if (fabs(changeDown)<1.0e30)
5319      sequenceIn_ = sequenceDown;
5320    alpha_ = alphaDown;
5321#ifdef CLP_DEBUG
5322    if ((handler_->logLevel()&32))
5323      printf("predicted way - dirout %d, change %g,%g theta %g\n",
5324             directionOut_,changeDown,changeUp,theta_);
5325#endif
5326  } else {
5327    theta_ = thetaUp;
5328    if (fabs(changeUp)<1.0e30)
5329      sequenceIn_ = sequenceUp;
5330    alpha_ = alphaUp;
5331    if (sequenceIn_!=sequenceOut_) {
5332#ifdef CLP_DEBUG
5333      if ((handler_->logLevel()&32))
5334        printf("opposite way - dirout %d, change %g,%g theta %g\n",
5335               directionOut_,changeDown,changeUp,theta_);
5336#endif
5337    } else {
5338#ifdef CLP_DEBUG
5339      if ((handler_->logLevel()&32))
5340        printf("opposite way to zero dj - dirout %d, change %g,%g theta %g\n",
5341               directionOut_,changeDown,changeUp,theta_);
5342#endif
5343    }
5344  }
5345  if (sequenceIn_>=0) {
5346    lowerIn_ = lower_[sequenceIn_];
5347    upperIn_ = upper_[sequenceIn_];
5348    valueIn_ = solution_[sequenceIn_];
5349    dualIn_ = dj_[sequenceIn_];
5350
5351    if (alpha_<0.0) {
5352      // as if from upper bound
5353      directionIn_=-1;
5354      upperIn_=valueIn_;
5355    } else {
5356      // as if from lower bound
5357      directionIn_=1;
5358      lowerIn_=valueIn_;
5359    }
5360  }
5361}
5362/*
5363   Row array has row part of pivot row
5364   Column array has column part.
5365   This is used in cleanup
5366*/
5367void
5368ClpSimplexDual::checkPossibleCleanup(CoinIndexedVector * rowArray,
5369                                        CoinIndexedVector * columnArray,
5370                                        double acceptablePivot)
5371{
5372  double * work;
5373  int number;
5374  int * which;
5375  int iSection;
5376
5377  double tolerance = dualTolerance_*1.001;
5378
5379  double thetaDown = 1.0e31;
5380  double thetaUp = 1.0e31;
5381  double bestAlphaDown = acceptablePivot*10.0;
5382  double bestAlphaUp = acceptablePivot*10.0;
5383  int sequenceDown =-1;
5384  int sequenceUp = -1;
5385
5386  double djSlack = dj_[pivotRow_];
5387  if (getRowStatus(pivotRow_)==basic)
5388    djSlack=COIN_DBL_MAX;
5389  if (fabs(djSlack)<tolerance)
5390    djSlack=0.0;
5391  int addSequence;
5392
5393  double alphaUp=0.0;
5394  double alphaDown=0.0;
5395  for (iSection=0;iSection<2;iSection++) {
5396
5397    int i;
5398    if (!iSection) {
5399      work = rowArray->denseVector();
5400      number = rowArray->getNumElements();
5401      which = rowArray->getIndices();
5402      addSequence = numberColumns_;
5403    } else {
5404      work = columnArray->denseVector();
5405      number = columnArray->getNumElements();
5406      which = columnArray->getIndices();
5407      addSequence = 0;
5408    }
5409   
5410    for (i=0;i<number;i++) {
5411      int iSequence = which[i];
5412      int iSequence2 = iSequence + addSequence;
5413      double alpha;
5414      double oldValue;
5415      double value;
5416
5417      switch(getStatus(iSequence2)) {
5418         
5419      case basic: 
5420        break;
5421      case ClpSimplex::isFixed:
5422        alpha = work[i];
5423        if (addSequence) {
5424          printf("possible - pivot row %d this %d\n",pivotRow_,iSequence);
5425          oldValue = dj_[iSequence2];
5426          if (alpha<=-acceptablePivot) {
5427            // might do other way
5428            value = oldValue+thetaUp*alpha;
5429            if (value<tolerance) {
5430              if (value<-tolerance||fabs(alpha)>bestAlphaUp) {
5431                thetaUp = -oldValue/alpha;
5432                bestAlphaUp = fabs(alpha);
5433                sequenceUp = iSequence2;
5434                alphaUp = alpha;
5435              }
5436            }
5437          } else if (alpha>=acceptablePivot) {
5438            // might do this way
5439            value = oldValue-thetaDown*alpha;
5440            if (value<tolerance) {
5441              if (value<-tolerance||fabs(alpha)>bestAlphaDown) {
5442                thetaDown = oldValue/alpha;
5443                bestAlphaDown = fabs(alpha);
5444                sequenceDown = iSequence2;
5445                alphaDown = alpha;
5446              }
5447            }
5448          }
5449        }
5450        break;
5451      case isFree:
5452      case superBasic:
5453        alpha = work[i];
5454        // dj must be effectively zero as dual feasible
5455        if (fabs(alpha)>bestAlphaUp) {
5456          thetaDown = 0.0;
5457          thetaUp = 0.0;
5458          bestAlphaDown = fabs(alpha);
5459          bestAlphaUp = bestAlphaDown;
5460          sequenceDown =iSequence2;
5461          sequenceUp = sequenceDown;
5462          alphaUp = alpha;
5463          alphaDown = alpha;
5464        }
5465        break;
5466      case atUpperBound:
5467        alpha = work[i];
5468        oldValue = dj_[iSequence2];
5469        if (alpha>=acceptablePivot) {
5470          // might do other way
5471          value = oldValue+thetaUp*alpha;
5472          if (value>-tolerance) {
5473            if (value>tolerance||fabs(alpha)>bestAlphaUp) {
5474              thetaUp = -oldValue/alpha;
5475              bestAlphaUp = fabs(alpha);
5476              sequenceUp = iSequence2;
5477              alphaUp = alpha;
5478            }
5479          }
5480        } else if (alpha<=-acceptablePivot) {
5481          // might do this way
5482          value = oldValue-thetaDown*alpha;
5483          if (value>-tolerance) {
5484            if (value>tolerance||fabs(alpha)>bestAlphaDown) {
5485              thetaDown = oldValue/alpha;
5486              bestAlphaDown = fabs(alpha);
5487              sequenceDown = iSequence2;
5488              alphaDown = alpha;
5489            }
5490          }
5491        }
5492        break;
5493      case atLowerBound:
5494        alpha = work[i];
5495        oldValue = dj_[iSequence2];
5496        if (alpha<=-acceptablePivot) {
5497          // might do other way
5498          value = oldValue+thetaUp*alpha;
5499          if (value<tolerance) {
5500            if (value<-tolerance||fabs(alpha)>bestAlphaUp) {
5501              thetaUp = -oldValue/alpha;
5502              bestAlphaUp = fabs(alpha);
5503              sequenceUp = iSequence2;
5504              alphaUp = alpha;
5505            }
5506          }
5507        } else if (alpha>=acceptablePivot) {
5508          // might do this way
5509          value = oldValue-thetaDown*alpha;
5510          if (value<tolerance) {
5511            if (value<-tolerance||fabs(alpha)>bestAlphaDown) {
5512              thetaDown = oldValue/alpha;
5513              bestAlphaDown = fabs(alpha);
5514              sequenceDown = iSequence2;
5515              alphaDown = alpha;
5516            }
5517          }
5518        }
5519        break;
5520      }
5521    }
5522  }
5523  thetaUp *= -1.0;
5524  // largest
5525  if (bestAlphaDown<bestAlphaUp) 
5526    sequenceDown =-1;
5527  else
5528    sequenceUp=-1;
5529
5530  sequenceIn_=-1;
5531 
5532  if (sequenceDown>=0) {
5533    theta_ = thetaDown;
5534    sequenceIn_ = sequenceDown;
5535    alpha_ = alphaDown;
5536#ifdef CLP_DEBUG
5537    if ((handler_->logLevel()&32))
5538      printf("predicted way - dirout %d, theta %g\n",
5539             directionOut_,theta_);
5540#endif
5541  } else if (sequenceUp>=0) {
5542    theta_ = thetaUp;
5543    sequenceIn_ = sequenceUp;
5544    alpha_ = alphaUp;
5545#ifdef CLP_DEBUG
5546    if ((handler_->logLevel()&32))
5547      printf("opposite way - dirout %d,theta %g\n",
5548             directionOut_,theta_);
5549#endif
5550  }
5551  if (sequenceIn_>=0) {
5552    lowerIn_ = lower_[sequenceIn_];
5553    upperIn_ = upper_[sequenceIn_];
5554    valueIn_ = solution_[sequenceIn_];
5555    dualIn_ = dj_[sequenceIn_];
5556
5557    if (alpha_<0.0) {
5558      // as if from upper bound
5559      directionIn_=-1;
5560      upperIn_=valueIn_;
5561    } else {
5562      // as if from lower bound
5563      directionIn_=1;
5564      lowerIn_=valueIn_;
5565    }
5566  }
5567}
5568/*
5569   This sees if we can move duals in dual values pass.
5570   This is done before any pivoting
5571*/
5572void ClpSimplexDual::doEasyOnesInValuesPass(double * dj)
5573{
5574  // Get column copy
5575  CoinPackedMatrix * columnCopy = matrix();
5576  // Get a row copy in standard format
5577  CoinPackedMatrix copy;
5578  copy.reverseOrderedCopyOf(*columnCopy);
5579  // get matrix data pointers
5580  const int * column = copy.getIndices();
5581  const CoinBigIndex * rowStart = copy.getVectorStarts();
5582  const int * rowLength = copy.getVectorLengths(); 
5583  const double * elementByRow = copy.getElements();
5584  double tolerance = dualTolerance_*1.001;
5585
5586  int iRow;
5587#ifdef CLP_DEBUG
5588  {
5589    double value5=0.0;
5590    int i;
5591    for (i=0;i<numberRows_+numberColumns_;i++) {
5592      if (dj[i]<-1.0e-6)
5593        value5 += dj[i]*upper_[i];
5594      else if (dj[i] >1.0e-6)
5595        value5 += dj[i]*lower_[i];
5596    }
5597    printf("Values objective Value before %g\n",value5);
5598  }
5599#endif
5600  // for scaled row
5601  double * scaled=NULL;
5602  if (rowScale_)
5603    scaled = new double[numberColumns_];
5604  for (iRow=0;iRow<numberRows_;iRow++) {
5605
5606    int iSequence = iRow + numberColumns_;
5607    double djBasic = dj[iSequence];
5608    if (getRowStatus(iRow)==basic&&fabs(djBasic)>tolerance) {
5609
5610      double changeUp ;
5611      // always -1 in pivot row
5612      if (djBasic>0.0) {
5613        // basic at lower bound
5614        changeUp = -lower_[iSequence];
5615      } else {
5616        // basic at upper bound
5617        changeUp = upper_[iSequence];
5618      }
5619      bool canMove=true;
5620      int i;
5621      const double * thisElements = elementByRow + rowStart[iRow]; 
5622      const int * thisIndices = column+rowStart[iRow];
5623      if (rowScale_) {
5624        assert (!auxiliaryModel_);
5625        // scale row
5626        double scale = rowScale_[iRow];
5627        for (i = 0;i<rowLength[iRow];i++) {
5628          int iColumn = thisIndices[i];
5629          double alpha = thisElements[i];
5630          scaled[i] = scale*alpha*columnScale_[iColumn];
5631        }
5632        thisElements = scaled;
5633      }
5634      for (i = 0;i<rowLength[iRow];i++) {
5635        int iColumn = thisIndices[i];
5636        double alpha = thisElements[i];
5637        double oldValue = dj[iColumn];;
5638        double value;
5639
5640        switch(getStatus(iColumn)) {
5641         
5642        case basic:
5643          if (dj[iColumn]<-tolerance&&
5644              fabs(solution_[iColumn]-upper_[iColumn])<1.0e-8) {
5645            // at ub
5646            changeUp += alpha*upper_[iColumn];
5647            // might do other way
5648            value = oldValue+djBasic*alpha;
5649            if (value>tolerance) 
5650              canMove=false;
5651          } else if (dj[iColumn]>tolerance&&
5652              fabs(solution_[iColumn]-lower_[iColumn])<1.0e-8) {
5653            changeUp += alpha*lower_[iColumn];
5654            // might do other way
5655            value = oldValue+djBasic*alpha;
5656            if (value<-tolerance)
5657              canMove=false;
5658          } else {
5659            canMove=false;
5660          }
5661          break;
5662        case ClpSimplex::isFixed:
5663          changeUp += alpha*upper_[iColumn];
5664          break;
5665        case isFree:
5666        case superBasic:
5667          canMove=false;
5668        break;
5669      case atUpperBound:
5670        changeUp += alpha*upper_[iColumn];
5671        // might do other way
5672        value = oldValue+djBasic*alpha;
5673        if (value>tolerance) 
5674          canMove=false;
5675        break;
5676      case atLowerBound:
5677        changeUp += alpha*lower_[iColumn];
5678        // might do other way
5679        value = oldValue+djBasic*alpha;
5680        if (value<-tolerance)
5681          canMove=false;
5682        break;
5683        }
5684      }
5685      if (canMove) {
5686        if (changeUp*djBasic>1.0e-12||fabs(changeUp)<1.0e-8) {
5687          // move
5688          for (i = 0;i<rowLength[iRow];i++) {
5689            int iColumn = thisIndices[i];
5690            double alpha = thisElements[i];
5691            dj[iColumn] += djBasic * alpha;
5692          }
5693          dj[iSequence]=0.0;
5694#ifdef CLP_DEBUG
5695          {
5696            double value5=0.0;
5697            int i;
5698            for (i=0;i<numberRows_+numberColumns_;i++) {
5699              if (dj[i]<-1.0e-6)
5700                value5 += dj[i]*upper_[i];
5701              else if (dj[i] >1.0e-6)
5702                value5 += dj[i]*lower_[i];
5703            }
5704            printf("Values objective Value after row %d old dj %g %g\n",
5705                   iRow,djBasic,value5);
5706          }
5707#endif
5708        }
5709      }
5710    }     
5711  }
5712  delete [] scaled;
5713}
5714int
5715ClpSimplexDual::nextSuperBasic()
5716{
5717  if (firstFree_>=0) {
5718    int returnValue=firstFree_;
5719    int iColumn=firstFree_+1;
5720    for (;iColumn<numberRows_+numberColumns_;iColumn++) {
5721      if (getStatus(iColumn)==isFree) 
5722        if (fabs(dj_[iColumn])>1.0e2*dualTolerance_) 
5723          break;
5724    }
5725    firstFree_ = iColumn;
5726    if (firstFree_==numberRows_+numberColumns_)
5727      firstFree_=-1;
5728    return returnValue;
5729  } else {
5730    return -1;
5731  }
5732}
5733void 
5734ClpSimplexDual::resetFakeBounds()
5735{
5736  double dummyChangeCost=0.0;
5737  changeBounds(true,rowArray_[2],dummyChangeCost);
5738  // throw away change
5739  for (int i=0;i<4;i++) 
5740    rowArray_[i]->clear();
5741}
Note: See TracBrowser for help on using the repository browser.