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

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

to pass miplib

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