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

Last change on this file since 1043 was 1043, checked in by forrest, 12 years ago

slightly more robust

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