source: stable/1.3/Clp/src/ClpSimplexDual.cpp @ 907

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

better cutoff checking

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 171.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,double & badFree)
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  badFree=0.0;
2418  for (int iSection=0;iSection<2;iSection++) {
2419
2420    int addSequence;
2421
2422    if (!iSection) {
2423      work = rowArray->denseVector();
2424      number = rowArray->getNumElements();
2425      which = rowArray->getIndices();
2426      reducedCost = rowReducedCost_;
2427      addSequence = numberColumns_;
2428    } else {
2429      work = columnArray->denseVector();
2430      number = columnArray->getNumElements();
2431      which = columnArray->getIndices();
2432      reducedCost = reducedCostWork_;
2433      addSequence = 0;
2434    }
2435   
2436    for (i=0;i<number;i++) {
2437      int iSequence = which[i];
2438      double alpha;
2439      double oldValue;
2440      double value;
2441      bool keep;
2442
2443      switch(getStatus(iSequence+addSequence)) {
2444         
2445      case basic:
2446      case ClpSimplex::isFixed:
2447        break;
2448      case isFree:
2449      case superBasic:
2450        alpha = work[i];
2451        bestPossible = CoinMax(bestPossible,fabs(alpha));
2452        oldValue = reducedCost[iSequence];
2453        // If free has to be very large - should come in via dualRow
2454        //if (getStatus(iSequence+addSequence)==isFree&&fabs(alpha)<1.0e-3)
2455        //break;
2456        if (oldValue>dualTolerance_) {
2457          keep = true;
2458        } else if (oldValue<-dualTolerance_) {
2459          keep = true;
2460        } else {
2461          if (fabs(alpha)>CoinMax(10.0*acceptablePivot,1.0e-5)) {
2462            keep = true;
2463          } else {
2464            keep=false;
2465            badFree=CoinMax(badFree,fabs(alpha));
2466          }
2467        }
2468        if (keep) {
2469          // free - choose largest
2470          if (fabs(alpha)>freePivot) {
2471            freePivot=fabs(alpha);
2472            sequenceIn_ = iSequence + addSequence;
2473            theta_=oldValue/alpha;
2474            alpha_=alpha;
2475          }
2476        }
2477        break;
2478      case atUpperBound:
2479        alpha = work[i];
2480        oldValue = reducedCost[iSequence];
2481        value = oldValue-tentativeTheta*alpha;
2482        //assert (oldValue<=dualTolerance_*1.0001);
2483        if (value>dualTolerance_) {
2484          bestPossible = CoinMax(bestPossible,-alpha);
2485          value = oldValue-upperTheta*alpha;
2486          if (value>dualTolerance_ && -alpha>=acceptablePivot) 
2487            upperTheta = (oldValue-dualTolerance_)/alpha;
2488          // add to list
2489          spare[numberRemaining]=alpha;
2490          index[numberRemaining++]=iSequence+addSequence;
2491        }
2492        break;
2493      case atLowerBound:
2494        alpha = work[i];
2495        oldValue = reducedCost[iSequence];
2496        value = oldValue-tentativeTheta*alpha;
2497        //assert (oldValue>=-dualTolerance_*1.0001);
2498        if (value<-dualTolerance_) {
2499          bestPossible = CoinMax(bestPossible,alpha);
2500          value = oldValue-upperTheta*alpha;
2501          if (value<-dualTolerance_ && alpha>=acceptablePivot) 
2502            upperTheta = (oldValue+dualTolerance_)/alpha;
2503          // add to list
2504          spare[numberRemaining]=alpha;
2505          index[numberRemaining++]=iSequence+addSequence;
2506        }
2507        break;
2508      }
2509    }
2510  }
2511  upperReturn = upperTheta;
2512  bestReturn = bestPossible;
2513  return numberRemaining;
2514}
2515/*
2516   Row array has row part of pivot row (as duals so sign may be switched)
2517   Column array has column part.
2518   This chooses pivot column.
2519   Spare array will be needed when we start getting clever.
2520   We will check for basic so spare array will never overflow.
2521   If necessary will modify costs
2522*/
2523double
2524ClpSimplexDual::dualColumn(CoinIndexedVector * rowArray,
2525                           CoinIndexedVector * columnArray,
2526                           CoinIndexedVector * spareArray,
2527                           CoinIndexedVector * spareArray2,
2528                           double acceptablePivot,
2529                           CoinBigIndex * dubiousWeights)
2530{
2531  int numberPossiblySwapped=0;
2532  int numberRemaining=0;
2533 
2534  double totalThru=0.0; // for when variables flip
2535  //double saveAcceptable=acceptablePivot;
2536  //acceptablePivot=1.0e-9;
2537
2538  double bestEverPivot=acceptablePivot;
2539  int lastSequence = -1;
2540  double lastPivot=0.0;
2541  double upperTheta;
2542  double newTolerance = dualTolerance_;
2543  //newTolerance = dualTolerance_+1.0e-6*dblParam_[ClpDualTolerance];
2544  // will we need to increase tolerance
2545  bool thisIncrease=false;
2546  // If we think we need to modify costs (not if something from broad sweep)
2547  bool modifyCosts=false;
2548  // Increase in objective due to swapping bounds (may be negative)
2549  double increaseInObjective=0.0;
2550
2551  // use spareArrays to put ones looked at in
2552  // we are going to flip flop between
2553  int iFlip = 0;
2554  // Possible list of pivots
2555  int interesting[2];
2556  // where possible swapped ones are
2557  int swapped[2];
2558  // for zeroing out arrays after
2559  int marker[2][2];
2560  // pivot elements
2561  double * array[2], * spare, * spare2;
2562  // indices
2563  int * indices[2], * index, * index2;
2564  spareArray2->clear();
2565  array[0] = spareArray->denseVector();
2566  indices[0] = spareArray->getIndices();
2567  spare = array[0];
2568  index = indices[0];
2569  array[1] = spareArray2->denseVector();
2570  indices[1] = spareArray2->getIndices();
2571  int i;
2572
2573  // initialize lists
2574  for (i=0;i<2;i++) {
2575    interesting[i]=0;
2576    swapped[i]=numberColumns_;
2577    marker[i][0]=0;
2578    marker[i][1]=numberColumns_;
2579  }
2580  /*
2581    First we get a list of possible pivots.  We can also see if the
2582    problem looks infeasible or whether we want to pivot in free variable.
2583    This may make objective go backwards but can only happen a finite
2584    number of times and I do want free variables basic.
2585
2586    Then we flip back and forth.  At the start of each iteration
2587    interesting[iFlip] should have possible candidates and swapped[iFlip]
2588    will have pivots if we decide to take a previous pivot.
2589    At end of each iteration interesting[1-iFlip] should have
2590    candidates if we go through this theta and swapped[1-iFlip]
2591    pivots if we don't go through.
2592
2593    At first we increase theta and see what happens.  We start
2594    theta at a reasonable guess.  If in right area then we do bit by bit.
2595
2596   */
2597
2598  // do first pass to get possibles
2599  upperTheta = 1.0e31;
2600  double bestPossible=0.0;
2601  double badFree=0.0;
2602  if (spareIntArray_[0]!=-1) {
2603    numberRemaining = dualColumn0(rowArray,columnArray,spareArray,
2604                                  acceptablePivot,upperTheta,bestPossible,badFree);
2605  } else {
2606    // already done
2607    numberRemaining = spareArray->getNumElements();
2608    spareArray->setNumElements(0);
2609    upperTheta = spareDoubleArray_[0];
2610    bestPossible = spareDoubleArray_[1];
2611    theta_ = spareDoubleArray_[2];
2612    alpha_ = spareDoubleArray_[3];
2613    sequenceIn_ = spareIntArray_[1];
2614  }
2615  // switch off
2616  spareIntArray_[0]=0;
2617  // We can also see if infeasible or pivoting on free
2618  double tentativeTheta = 1.0e25;
2619  interesting[0]=numberRemaining;
2620  marker[0][0] = numberRemaining;
2621
2622  if (!numberRemaining&&sequenceIn_<0)
2623    return 0.0; // Looks infeasible
2624
2625  // If sum of bad small pivots too much
2626#define MORE_CAREFUL
2627#ifdef MORE_CAREFUL
2628  bool badSumPivots=false;
2629#endif
2630  if (sequenceIn_>=0) {
2631    // free variable - always choose
2632  } else {
2633
2634    theta_=1.0e50;
2635    // now flip flop between spare arrays until reasonable theta
2636    tentativeTheta = CoinMax(10.0*upperTheta,1.0e-7);
2637   
2638    // loops increasing tentative theta until can't go through
2639   
2640    while (tentativeTheta < 1.0e22) {
2641      double thruThis = 0.0;
2642     
2643      double bestPivot=acceptablePivot;
2644      int bestSequence=-1;
2645     
2646      numberPossiblySwapped = numberColumns_;
2647      numberRemaining = 0;
2648
2649      upperTheta = 1.0e50;
2650
2651      spare = array[iFlip];
2652      index = indices[iFlip];
2653      spare2 = array[1-iFlip];
2654      index2 = indices[1-iFlip];
2655
2656      // try 3 different ways
2657      // 1 bias increase by ones with slightly wrong djs
2658      // 2 bias by all
2659      // 3 bias by all - tolerance
2660#define TRYBIAS 3
2661
2662
2663      double increaseInThis=0.0; //objective increase in this loop
2664     
2665      for (i=0;i<interesting[iFlip];i++) {
2666        int iSequence = index[i];
2667        double alpha = spare[i];
2668        double oldValue = dj_[iSequence];
2669        double value = oldValue-tentativeTheta*alpha;
2670
2671        if (alpha < 0.0) {
2672          //at upper bound
2673          if (value>newTolerance) {
2674            double range = upper_[iSequence] - lower_[iSequence];
2675            thruThis -= range*alpha;
2676#if TRYBIAS==1
2677            if (oldValue>0.0)
2678              increaseInThis -= oldValue*range;
2679#elif TRYBIAS==2
2680            increaseInThis -= oldValue*range;
2681#else
2682            increaseInThis -= (oldValue+dualTolerance_)*range;
2683#endif
2684            // goes on swapped list (also means candidates if too many)
2685            spare2[--numberPossiblySwapped]=alpha;
2686            index2[numberPossiblySwapped]=iSequence;
2687            if (fabs(alpha)>bestPivot) {
2688              bestPivot=fabs(alpha);
2689              bestSequence=numberPossiblySwapped;
2690            }
2691          } else {
2692            value = oldValue-upperTheta*alpha;
2693            if (value>newTolerance && -alpha>=acceptablePivot) 
2694              upperTheta = (oldValue-newTolerance)/alpha;
2695            spare2[numberRemaining]=alpha;
2696            index2[numberRemaining++]=iSequence;
2697          }
2698        } else {
2699          // at lower bound
2700          if (value<-newTolerance) {
2701            double range = upper_[iSequence] - lower_[iSequence];
2702            thruThis += range*alpha;
2703            //?? is this correct - and should we look at good ones
2704#if TRYBIAS==1
2705            if (oldValue<0.0)
2706              increaseInThis += oldValue*range;
2707#elif TRYBIAS==2
2708            increaseInThis += oldValue*range;
2709#else
2710            increaseInThis += (oldValue-dualTolerance_)*range;
2711#endif
2712            // goes on swapped list (also means candidates if too many)
2713            spare2[--numberPossiblySwapped]=alpha;
2714            index2[numberPossiblySwapped]=iSequence;
2715            if (fabs(alpha)>bestPivot) {
2716              bestPivot=fabs(alpha);
2717              bestSequence=numberPossiblySwapped;
2718            }
2719          } else {
2720            value = oldValue-upperTheta*alpha;
2721            if (value<-newTolerance && alpha>=acceptablePivot) 
2722              upperTheta = (oldValue+newTolerance)/alpha;
2723            spare2[numberRemaining]=alpha;
2724            index2[numberRemaining++]=iSequence;
2725          }
2726        }
2727      }
2728      swapped[1-iFlip]=numberPossiblySwapped;
2729      interesting[1-iFlip]=numberRemaining;
2730      marker[1-iFlip][0]= CoinMax(marker[1-iFlip][0],numberRemaining);
2731      marker[1-iFlip][1]= CoinMin(marker[1-iFlip][1],numberPossiblySwapped);
2732     
2733      if (totalThru+thruThis>=fabs(dualOut_)||
2734          increaseInObjective+increaseInThis<0.0) {
2735        // We should be pivoting in this batch
2736        // so compress down to this lot
2737        numberRemaining=0;
2738        for (i=numberColumns_-1;i>=swapped[1-iFlip];i--) {
2739          spare[numberRemaining]=spare2[i];
2740          index[numberRemaining++]=index2[i];
2741        }
2742        interesting[iFlip]=numberRemaining;
2743        int iTry;
2744#define MAXTRY 100
2745        // first get ratio with tolerance
2746        for (iTry=0;iTry<MAXTRY;iTry++) {
2747         
2748          upperTheta=1.0e50;
2749          numberPossiblySwapped = numberColumns_;
2750          numberRemaining = 0;
2751
2752          increaseInThis=0.0; //objective increase in this loop
2753
2754          thruThis=0.0;
2755
2756          spare = array[iFlip];
2757          index = indices[iFlip];
2758          spare2 = array[1-iFlip];
2759          index2 = indices[1-iFlip];
2760          for (i=0;i<interesting[iFlip];i++) {
2761            int iSequence=index[i];
2762            double alpha=spare[i];
2763            double oldValue = dj_[iSequence];
2764            double value = oldValue-upperTheta*alpha;
2765           
2766            if (alpha < 0.0) {
2767              //at upper bound
2768              if (value>newTolerance) {
2769                if (-alpha>=acceptablePivot) {
2770                  upperTheta = (oldValue-newTolerance)/alpha;
2771                }
2772              }
2773            } else {
2774              // at lower bound
2775              if (value<-newTolerance) {
2776                if (alpha>=acceptablePivot) {
2777                  upperTheta = (oldValue+newTolerance)/alpha;
2778                }
2779              }
2780            }
2781          }
2782          bestPivot=acceptablePivot;
2783          sequenceIn_=-1;
2784#ifdef DUBIOUS_WEIGHTS
2785          double bestWeight=COIN_DBL_MAX;
2786#endif
2787          double largestPivot=acceptablePivot;
2788          // now choose largest and sum all ones which will go through
2789          //printf("XX it %d number %d\n",numberIterations_,interesting[iFlip]);
2790  // Sum of bad small pivots
2791#ifdef MORE_CAREFUL
2792          double sumBadPivots=0.0;
2793          badSumPivots=false;
2794#endif
2795          // Make sure upperTheta will work (-O2 and above gives problems)
2796          upperTheta *= 1.0000000001;
2797          for (i=0;i<interesting[iFlip];i++) {
2798            int iSequence=index[i];
2799            double alpha=spare[i];
2800            double value = dj_[iSequence]-upperTheta*alpha;
2801            double badDj=0.0;
2802           
2803            bool addToSwapped=false;
2804           
2805            if (alpha < 0.0) {
2806              //at upper bound
2807              if (value>=0.0) { 
2808                addToSwapped=true;
2809#if TRYBIAS==1
2810                badDj = -CoinMax(dj_[iSequence],0.0);
2811#elif TRYBIAS==2
2812                badDj = -dj_[iSequence];
2813#else
2814                badDj = -dj_[iSequence]-dualTolerance_;
2815#endif
2816              }
2817            } else {
2818              // at lower bound
2819              if (value<=0.0) {
2820                addToSwapped=true;
2821#if TRYBIAS==1
2822                badDj = CoinMin(dj_[iSequence],0.0);
2823#elif TRYBIAS==2
2824                badDj = dj_[iSequence];
2825#else
2826                badDj = dj_[iSequence]-dualTolerance_;
2827#endif
2828              }
2829            }
2830            if (!addToSwapped) {
2831              // add to list of remaining
2832              spare2[numberRemaining]=alpha;
2833              index2[numberRemaining++]=iSequence;
2834            } else {
2835              // add to list of swapped
2836              spare2[--numberPossiblySwapped]=alpha;
2837              index2[numberPossiblySwapped]=iSequence;
2838              // select if largest pivot
2839              bool take=false;
2840              double absAlpha = fabs(alpha);
2841#ifdef DUBIOUS_WEIGHTS
2842              // User could do anything to break ties here
2843              double weight;
2844              if (dubiousWeights)
2845                weight=dubiousWeights[iSequence];
2846              else
2847                weight=1.0;
2848              weight += CoinDrand48()*1.0e-2;
2849              if (absAlpha>2.0*bestPivot) {
2850                take=true;
2851              } else if (absAlpha>largestPivot) {
2852                // could multiply absAlpha and weight
2853                if (weight*bestPivot<bestWeight*absAlpha)
2854                  take=true;
2855              }
2856#else
2857              if (absAlpha>bestPivot) 
2858                take=true;
2859#endif
2860#ifdef MORE_CAREFUL
2861              if (absAlpha<acceptablePivot&&upperTheta<1.0e20) {
2862                if (alpha < 0.0) {
2863                  //at upper bound
2864                  if (value>dualTolerance_) {
2865                    double gap=upper_[iSequence]-lower_[iSequence];
2866                    if (gap<1.0e20)
2867                      sumBadPivots += value*gap; 
2868                    else
2869                      sumBadPivots += 1.0e20;
2870                    //printf("bad %d alpha %g dj at upper %g\n",
2871                    //     iSequence,alpha,value);
2872                  }
2873                } else {
2874                  //at lower bound
2875                  if (value<-dualTolerance_) {
2876                    double gap=upper_[iSequence]-lower_[iSequence];
2877                    if (gap<1.0e20)
2878                      sumBadPivots -= value*gap; 
2879                    else
2880                      sumBadPivots += 1.0e20;
2881                    //printf("bad %d alpha %g dj at lower %g\n",
2882                    //     iSequence,alpha,value);
2883                  }
2884                }
2885              }
2886#endif
2887#ifdef FORCE_FOLLOW
2888              if (iSequence==force_in) {
2889                printf("taking %d - alpha %g best %g\n",force_in,absAlpha,largestPivot);
2890                take=true;
2891              }
2892#endif
2893              if (take) {
2894                sequenceIn_ = numberPossiblySwapped;
2895                bestPivot =  absAlpha;
2896                theta_ = dj_[iSequence]/alpha;
2897                largestPivot = CoinMax(largestPivot,0.5*bestPivot);
2898#ifdef DUBIOUS_WEIGHTS
2899                bestWeight = weight;
2900#endif
2901                //printf(" taken seq %d alpha %g weight %d\n",
2902                //   iSequence,absAlpha,dubiousWeights[iSequence]);
2903              } else {
2904                //printf(" not taken seq %d alpha %g weight %d\n",
2905                //   iSequence,absAlpha,dubiousWeights[iSequence]);
2906              }
2907              double range = upper_[iSequence] - lower_[iSequence];
2908              thruThis += range*fabs(alpha);
2909              increaseInThis += badDj*range;
2910            }
2911          }
2912#ifdef MORE_CAREFUL
2913          // If we have done pivots and things look bad set alpha_ 0.0 to force factorization
2914          if (sumBadPivots>1.0e4) {
2915            if (handler_->logLevel()>1)
2916            printf("maybe forcing re-factorization - sum %g  %d pivots\n",sumBadPivots,
2917                   factorization_->pivots());
2918            badSumPivots=true;
2919          }
2920#endif
2921          swapped[1-iFlip]=numberPossiblySwapped;
2922          interesting[1-iFlip]=numberRemaining;
2923          marker[1-iFlip][0]= CoinMax(marker[1-iFlip][0],numberRemaining);
2924          marker[1-iFlip][1]= CoinMin(marker[1-iFlip][1],numberPossiblySwapped);
2925          // If we stop now this will be increase in objective (I think)
2926          double increase = (fabs(dualOut_)-totalThru)*theta_;
2927          increase += increaseInObjective;
2928          if (theta_<0.0)
2929            thruThis += fabs(dualOut_); // force using this one
2930          if (increaseInObjective<0.0&&increase<0.0&&lastSequence>=0) {
2931            // back
2932            // We may need to be more careful - we could do by
2933            // switch so we always do fine grained?
2934            bestPivot=0.0;
2935          } else {
2936            // add in
2937            totalThru += thruThis;
2938            increaseInObjective += increaseInThis;
2939          }
2940          if (bestPivot<0.1*bestEverPivot&&
2941              bestEverPivot>1.0e-6&&
2942              (bestPivot<1.0e-3||totalThru*2.0>fabs(dualOut_))) {
2943            // back to previous one
2944            sequenceIn_=lastSequence;
2945            // swap regions
2946            iFlip = 1-iFlip;
2947            break;
2948          } else if (sequenceIn_==-1&&upperTheta>largeValue_) {
2949            if (lastPivot>acceptablePivot) {
2950              // back to previous one
2951              sequenceIn_=lastSequence;
2952              // swap regions
2953              iFlip = 1-iFlip;
2954            } else {
2955              // can only get here if all pivots too small
2956            }
2957            break;
2958          } else if (totalThru>=fabs(dualOut_)) {
2959            modifyCosts=true; // fine grain - we can modify costs
2960            break; // no point trying another loop
2961          } else {
2962            lastSequence=sequenceIn_;
2963            if (bestPivot>bestEverPivot)
2964              bestEverPivot=bestPivot;
2965            iFlip = 1 -iFlip;
2966            modifyCosts=true; // fine grain - we can modify costs
2967          }
2968        }
2969        if (iTry==MAXTRY)
2970          iFlip = 1-iFlip; // flip back
2971        break;
2972      } else {
2973        // skip this lot
2974        if (bestPivot>1.0e-3||bestPivot>bestEverPivot) {
2975          bestEverPivot=bestPivot;
2976          lastSequence=bestSequence;
2977        } else {
2978          // keep old swapped
2979          CoinMemcpyN(array[iFlip]+swapped[iFlip],
2980                 numberColumns_-swapped[iFlip],array[1-iFlip]+swapped[iFlip]);
2981          CoinMemcpyN(indices[iFlip]+swapped[iFlip],
2982                 numberColumns_-swapped[iFlip],indices[1-iFlip]+swapped[iFlip]);
2983          marker[1-iFlip][1] = CoinMin(marker[1-iFlip][1],swapped[iFlip]);
2984          swapped[1-iFlip]=swapped[iFlip];
2985        }
2986        increaseInObjective += increaseInThis;
2987        iFlip = 1 - iFlip; // swap regions
2988        tentativeTheta = 2.0*upperTheta;
2989        totalThru += thruThis;
2990      }
2991    }
2992   
2993    // can get here without sequenceIn_ set but with lastSequence
2994    if (sequenceIn_<0&&lastSequence>=0) {
2995      // back to previous one
2996      sequenceIn_=lastSequence;
2997      // swap regions
2998      iFlip = 1-iFlip;
2999    }
3000   
3001#define MINIMUMTHETA 1.0e-18
3002    // Movement should be minimum for anti-degeneracy - unless
3003    // fixed variable out
3004    double minimumTheta;
3005    if (upperOut_>lowerOut_)
3006      minimumTheta=MINIMUMTHETA;
3007    else
3008      minimumTheta=0.0;
3009    if (sequenceIn_>=0) {
3010      // at this stage sequenceIn_ is just pointer into index array
3011      // flip just so we can use iFlip
3012      iFlip = 1 -iFlip;
3013      spare = array[iFlip];
3014      index = indices[iFlip];
3015      double oldValue;
3016      alpha_ = spare[sequenceIn_];
3017      sequenceIn_ = indices[iFlip][sequenceIn_];
3018      oldValue = dj_[sequenceIn_];
3019      theta_ = CoinMax(oldValue/alpha_,0.0);
3020      if (theta_<minimumTheta&&fabs(alpha_)<1.0e5&&1) {
3021        // can't pivot to zero
3022#if 0
3023        if (oldValue-minimumTheta*alpha_>=-dualTolerance_) {
3024          theta_=minimumTheta;
3025        } else if (oldValue-minimumTheta*alpha_>=-newTolerance) {
3026          theta_=minimumTheta;
3027          thisIncrease=true;
3028        } else {
3029          theta_=CoinMax((oldValue+newTolerance)/alpha_,0.0);
3030          thisIncrease=true;
3031        }
3032#else
3033        theta_=minimumTheta;
3034#endif
3035      }
3036      // may need to adjust costs so all dual feasible AND pivoted is exactly 0
3037      //int costOffset = numberRows_+numberColumns_;
3038      if (modifyCosts) {
3039        int i;
3040        for (i=numberColumns_-1;i>=swapped[iFlip];i--) {
3041          int iSequence=index[i];
3042          double alpha=spare[i];
3043          double value = dj_[iSequence]-theta_*alpha;
3044           
3045          // can't be free here
3046         
3047          if (alpha < 0.0) {
3048            //at upper bound
3049            if (value>dualTolerance_) {
3050              thisIncrease=true;
3051#define MODIFYCOST 2
3052#if MODIFYCOST
3053              // modify cost to hit new tolerance
3054              double modification = alpha*theta_-dj_[iSequence]
3055                +newTolerance;
3056              if ((specialOptions_&(2048+4096+16384))!=0) {
3057                if ((specialOptions_&16384)!=0) {
3058                  if (fabs(modification)<1.0e-8)
3059                    modification=0.0;
3060                } else if ((specialOptions_&2048)!=0) {
3061                  if (fabs(modification)<1.0e-10)
3062                    modification=0.0;
3063                } else {
3064                  if (fabs(modification)<1.0e-12)
3065                    modification=0.0;
3066                }
3067              }
3068              dj_[iSequence] += modification;
3069              cost_[iSequence] +=  modification;
3070              if (modification)
3071                numberChanged_ ++; // Say changed costs
3072              //cost_[iSequence+costOffset] += modification; // save change
3073#endif
3074            }
3075          } else {
3076            // at lower bound
3077            if (-value>dualTolerance_) {
3078              thisIncrease=true;
3079#if MODIFYCOST
3080              // modify cost to hit new tolerance
3081              double modification = alpha*theta_-dj_[iSequence]
3082                -newTolerance;
3083              //modification = CoinMax(modification,-dualTolerance_);
3084              //assert (fabs(modification)<1.0e-7);
3085              if ((specialOptions_&(2048+4096))!=0) {
3086                if ((specialOptions_&2048)!=0) {
3087                  if (fabs(modification)<1.0e-10)
3088                    modification=0.0;
3089                } else {
3090                  if (fabs(modification)<1.0e-12)
3091                    modification=0.0;
3092                }
3093              }
3094              dj_[iSequence] += modification;
3095              cost_[iSequence] +=  modification;
3096              if (modification)
3097                numberChanged_ ++; // Say changed costs
3098              //cost_[iSequence+costOffset] += modification; // save change
3099#endif
3100            }
3101          }
3102        }
3103      }
3104    }
3105  }
3106
3107  if (sequenceIn_>=0) {
3108#ifdef MORE_CAREFUL
3109    // If we have done pivots and things look bad set alpha_ 0.0 to force factorization
3110    if (badSumPivots&&factorization_->pivots()) {
3111      if (handler_->logLevel()>1)
3112        printf("forcing re-factorization\n");
3113      alpha_=0.0;
3114    }
3115    if (fabs(theta_*badFree)>10.0*dualTolerance_&&factorization_->pivots()) {
3116      if (handler_->logLevel()>1)
3117        printf("forcing re-factorizationon free\n");
3118      alpha_=0.0;
3119    }
3120#endif
3121    lowerIn_ = lower_[sequenceIn_];
3122    upperIn_ = upper_[sequenceIn_];
3123    valueIn_ = solution_[sequenceIn_];
3124    dualIn_ = dj_[sequenceIn_];
3125
3126    if (numberTimesOptimal_) {
3127      // can we adjust cost back closer to original
3128      //*** add coding
3129    }
3130#if MODIFYCOST>1   
3131    // modify cost to hit zero exactly
3132    // so (dualIn_+modification)==theta_*alpha_
3133    double modification = theta_*alpha_-dualIn_;
3134    if ((specialOptions_&(2048+4096))!=0) {
3135      if ((specialOptions_&16384)!=0) {
3136        // in fast dual
3137        if (fabs(modification)<1.0e-7)
3138          modification=0.0;
3139      } else if ((specialOptions_&2048)!=0) {
3140        if (fabs(modification)<1.0e-10)
3141          modification=0.0;
3142      } else {
3143        if (fabs(modification)<1.0e-12)
3144          modification=0.0;
3145      }
3146    }
3147    dualIn_ += modification;
3148    dj_[sequenceIn_]=dualIn_;
3149    cost_[sequenceIn_] += modification;
3150    if (modification)
3151      numberChanged_ ++; // Say changed costs
3152    //int costOffset = numberRows_+numberColumns_;
3153    //cost_[sequenceIn_+costOffset] += modification; // save change
3154    //assert (fabs(modification)<1.0e-6);
3155#ifdef CLP_DEBUG
3156    if ((handler_->logLevel()&32)&&fabs(modification)>1.0e-15)
3157      printf("exact %d new cost %g, change %g\n",sequenceIn_,
3158             cost_[sequenceIn_],modification);
3159#endif
3160#endif
3161   
3162    if (alpha_<0.0) {
3163      // as if from upper bound
3164      directionIn_=-1;
3165      upperIn_=valueIn_;
3166    } else {
3167      // as if from lower bound
3168      directionIn_=1;
3169      lowerIn_=valueIn_;
3170    }
3171  }
3172  //if (thisIncrease)
3173  //dualTolerance_+= 1.0e-6*dblParam_[ClpDualTolerance];
3174
3175  // clear arrays
3176
3177  for (i=0;i<2;i++) {
3178    CoinZeroN(array[i],marker[i][0]);
3179    CoinZeroN(array[i]+marker[i][1],numberColumns_-marker[i][1]);
3180  }
3181  return bestPossible;
3182}
3183#ifdef CLP_ALL_ONE_FILE
3184#undef MAXTRY
3185#endif
3186/* Checks if tentative optimal actually means unbounded
3187   Returns -3 if not, 2 if is unbounded */
3188int 
3189ClpSimplexDual::checkUnbounded(CoinIndexedVector * ray,
3190                                   CoinIndexedVector * spare,
3191                                   double changeCost)
3192{
3193  int status=2; // say unbounded
3194  factorization_->updateColumn(spare,ray);
3195  // get reduced cost
3196  int i;
3197  int number=ray->getNumElements();
3198  int * index = ray->getIndices();
3199  double * array = ray->denseVector();
3200  for (i=0;i<number;i++) {
3201    int iRow=index[i];
3202    int iPivot=pivotVariable_[iRow];
3203    changeCost -= cost(iPivot)*array[iRow];
3204  }
3205  double way;
3206  if (changeCost>0.0) {
3207    //try going down
3208    way=1.0;
3209  } else if (changeCost<0.0) {
3210    //try going up
3211    way=-1.0;
3212  } else {
3213#ifdef CLP_DEBUG
3214    printf("can't decide on up or down\n");
3215#endif
3216    way=0.0;
3217    status=-3;
3218  }
3219  double movement=1.0e10*way; // some largish number
3220  double zeroTolerance = 1.0e-14*dualBound_;
3221  for (i=0;i<number;i++) {
3222    int iRow=index[i];
3223    int iPivot=pivotVariable_[iRow];
3224    double arrayValue = array[iRow];
3225    if (fabs(arrayValue)<zeroTolerance)
3226      arrayValue=0.0;
3227    double newValue=solution(iPivot)+movement*arrayValue;
3228    if (newValue>upper(iPivot)+primalTolerance_||
3229        newValue<lower(iPivot)-primalTolerance_)
3230      status=-3; // not unbounded
3231  }
3232  if (status==2) {
3233    // create ray
3234    delete [] ray_;
3235    ray_ = new double [numberColumns_];
3236    CoinZeroN(ray_,numberColumns_);
3237    for (i=0;i<number;i++) {
3238      int iRow=index[i];
3239      int iPivot=pivotVariable_[iRow];
3240      double arrayValue = array[iRow];
3241      if (iPivot<numberColumns_&&fabs(arrayValue)>=zeroTolerance)
3242        ray_[iPivot] = way* array[iRow];
3243    }
3244  }
3245  ray->clear();
3246  return status;
3247}
3248/* Checks if finished.  Updates status */
3249void 
3250ClpSimplexDual::statusOfProblemInDual(int & lastCleaned,int type,
3251                                      double * givenDuals, ClpDataSave & saveData,
3252                                      int ifValuesPass)
3253{
3254  // If lots of iterations then adjust costs if large ones
3255  if (numberIterations_>4*(numberRows_+numberColumns_)&&objectiveScale_==1.0) {
3256    double largest=0.0;
3257    for (int i=0;i<numberRows_;i++) {
3258      int iColumn = pivotVariable_[i];
3259      largest = CoinMax(largest,fabs(cost_[iColumn]));
3260    }
3261    if (largest>1.0e6) {
3262      objectiveScale_ = 1.0e6/largest;
3263      for (int i=0;i<numberRows_+numberColumns_;i++)
3264        cost_[i] *= objectiveScale_;
3265    }
3266  }
3267  bool normalType=true;
3268  int numberPivots = factorization_->pivots();
3269  double realDualInfeasibilities=0.0;
3270  if (type==2) {
3271    // trouble - restore solution
3272    CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
3273    CoinMemcpyN(savedSolution_+numberColumns_ ,
3274           numberRows_,rowActivityWork_);
3275    CoinMemcpyN(savedSolution_ ,
3276           numberColumns_,columnActivityWork_);
3277    // restore extra stuff
3278    int dummy;
3279    matrix_->generalExpanded(this,6,dummy);
3280    forceFactorization_=1; // a bit drastic but ..
3281    changeMade_++; // say something changed
3282  }
3283  int tentativeStatus = problemStatus_;
3284  double changeCost;
3285  bool unflagVariables = true;
3286  if (problemStatus_>-3||factorization_->pivots()) {
3287    // factorize
3288    // later on we will need to recover from singularities
3289    // also we could skip if first time
3290    // save dual weights
3291    dualRowPivot_->saveWeights(this,1);
3292    if (type) {
3293      // is factorization okay?
3294      if (internalFactorize(1)) {
3295        // no - restore previous basis
3296        unflagVariables = false;
3297        assert (type==1);
3298        changeMade_++; // say something changed
3299        // Keep any flagged variables
3300        int i;
3301        for (i=0;i<numberRows_+numberColumns_;i++) {
3302          if (flagged(i))
3303            saveStatus_[i] |= 64; //say flagged
3304        }
3305        CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
3306        CoinMemcpyN(savedSolution_+numberColumns_ ,
3307               numberRows_,rowActivityWork_);
3308        CoinMemcpyN(savedSolution_ ,
3309               numberColumns_,columnActivityWork_);
3310        // restore extra stuff
3311        int dummy;
3312        matrix_->generalExpanded(this,6,dummy);
3313        // get correct bounds on all variables
3314        double dummyChangeCost=0.0;
3315        changeBounds(true,rowArray_[2],dummyChangeCost);
3316        // throw away change
3317        for (i=0;i<4;i++) 
3318          rowArray_[i]->clear();
3319        // need to reject something
3320        char x = isColumn(sequenceOut_) ? 'C' :'R';
3321        handler_->message(CLP_SIMPLEX_FLAG,messages_)
3322          <<x<<sequenceWithin(sequenceOut_)
3323          <<CoinMessageEol;
3324        setFlagged(sequenceOut_);
3325        progress_->clearBadTimes();
3326       
3327        // Go to safe
3328        factorization_->pivotTolerance(0.99);
3329        forceFactorization_=1; // a bit drastic but ..
3330        type = 2;
3331        //assert (internalFactorize(1)==0);
3332        if (internalFactorize(1)) {
3333          CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
3334          CoinMemcpyN(savedSolution_+numberColumns_ ,
3335                 numberRows_,rowActivityWork_);
3336          CoinMemcpyN(savedSolution_ ,
3337                 numberColumns_,columnActivityWork_);
3338          // restore extra stuff
3339          int dummy;
3340          matrix_->generalExpanded(this,6,dummy);
3341          // debug
3342          int returnCode = internalFactorize(1);
3343          while (returnCode) {
3344            // ouch
3345            // switch off dense
3346            int saveDense = factorization_->denseThreshold();
3347            factorization_->setDenseThreshold(0);
3348            // Go to safe
3349            factorization_->pivotTolerance(0.99);
3350            // make sure will do safe factorization
3351            pivotVariable_[0]=-1;
3352            returnCode=internalFactorize(2);
3353            factorization_->setDenseThreshold(saveDense);
3354          }
3355        }
3356      }
3357    }
3358    if (problemStatus_!=-4||factorization_->pivots()>10)
3359      problemStatus_=-3;
3360  }
3361  // at this stage status is -3 or -4 if looks infeasible
3362  // get primal and dual solutions
3363  gutsOfSolution(givenDuals,NULL);
3364  // If bad accuracy treat as singular
3365  if ((largestPrimalError_>1.0e15||largestDualError_>1.0e15)&&numberIterations_) {
3366    // restore previous basis
3367    unflagVariables = false;
3368    changeMade_++; // say something changed
3369    // Keep any flagged variables
3370    int i;
3371    for (i=0;i<numberRows_+numberColumns_;i++) {
3372      if (flagged(i))
3373        saveStatus_[i] |= 64; //say flagged
3374    }
3375    CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
3376    CoinMemcpyN(savedSolution_+numberColumns_ ,
3377           numberRows_,rowActivityWork_);
3378    CoinMemcpyN(savedSolution_ ,
3379           numberColumns_,columnActivityWork_);
3380    // restore extra stuff
3381    int dummy;
3382    matrix_->generalExpanded(this,6,dummy);
3383    // get correct bounds on all variables
3384    //double dummyChangeCost=0.0;
3385    //changeBounds(true,rowArray_[2],dummyChangeCost);
3386    // throw away change
3387    for (i=0;i<4;i++) 
3388      rowArray_[i]->clear();
3389    // need to reject something
3390    char x = isColumn(sequenceOut_) ? 'C' :'R';
3391    handler_->message(CLP_SIMPLEX_FLAG,messages_)
3392      <<x<<sequenceWithin(sequenceOut_)
3393      <<CoinMessageEol;
3394    setFlagged(sequenceOut_);
3395    progress_->clearBadTimes();
3396   
3397    // Go to safer
3398    double newTolerance = CoinMin(1.1*factorization_->pivotTolerance(),0.99);
3399    factorization_->pivotTolerance(newTolerance);
3400    forceFactorization_=1; // a bit drastic but ..
3401    type = 2;
3402    //assert (internalFactorize(1)==0);
3403    if (internalFactorize(1)) {
3404      CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
3405      CoinMemcpyN(savedSolution_+numberColumns_ ,
3406              numberRows_,rowActivityWork_);
3407      CoinMemcpyN(savedSolution_ ,
3408           numberColumns_,columnActivityWork_);
3409      // restore extra stuff
3410      int dummy;
3411      matrix_->generalExpanded(this,6,dummy);
3412      // debug
3413      int returnCode = internalFactorize(1);
3414      while (returnCode) {
3415        // ouch
3416        // switch off dense
3417        int saveDense = factorization_->denseThreshold();
3418        factorization_->setDenseThreshold(0);
3419        // Go to safe
3420        factorization_->pivotTolerance(0.99);
3421        // make sure will do safe factorization
3422        pivotVariable_[0]=-1;
3423        returnCode=internalFactorize(2);
3424        factorization_->setDenseThreshold(saveDense);
3425      }
3426    }
3427    // get primal and dual solutions
3428    gutsOfSolution(givenDuals,NULL);
3429  } else if (largestPrimalError_<1.0e-7&&largestDualError_<1.0e-7) {
3430    // Can reduce tolerance
3431    double newTolerance = CoinMax(0.99*factorization_->pivotTolerance(),saveData.pivotTolerance_);
3432    factorization_->pivotTolerance(newTolerance);
3433  } 
3434  // Double check infeasibility if no action
3435  if (progress_->lastIterationNumber(0)==numberIterations_) {
3436    if (dualRowPivot_->looksOptimal()) {
3437      numberPrimalInfeasibilities_ = 0;
3438      sumPrimalInfeasibilities_ = 0.0;
3439    }
3440#if 1
3441  } else {
3442    double thisObj = objectiveValue_;
3443    double lastObj = progress_->lastObjective(0);
3444    if (lastObj>thisObj+1.0e-3*CoinMax(fabs(thisObj),fabs(lastObj))+1.0
3445        &&!ifValuesPass&&firstFree_<0) {
3446      int maxFactor = factorization_->maximumPivots();
3447      if (maxFactor>10&&numberPivots>1) {
3448        //printf("lastobj %g thisobj %g\n",lastObj,thisObj);
3449        //if (forceFactorization_<0)
3450        //forceFactorization_= maxFactor;
3451        //forceFactorization_ = CoinMax(1,(forceFactorization_>>1));
3452        forceFactorization_=1;
3453        //printf("Reducing factorization frequency - bad backwards\n");
3454        unflagVariables = false;
3455        changeMade_++; // say something changed
3456        CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
3457        CoinMemcpyN(savedSolution_+numberColumns_ ,
3458                numberRows_,rowActivityWork_);
3459        CoinMemcpyN(savedSolution_ ,
3460                numberColumns_,columnActivityWork_);
3461        // restore extra stuff
3462        int dummy;
3463        matrix_->generalExpanded(this,6,dummy);
3464        // get correct bounds on all variables
3465        double dummyChangeCost=0.0;
3466        changeBounds(true,rowArray_[2],dummyChangeCost);
3467        // throw away change
3468        for (int i=0;i<4;i++) 
3469          rowArray_[i]->clear();
3470        if(factorization_->pivotTolerance()<0.2)
3471          factorization_->pivotTolerance(0.2);
3472        if (internalFactorize(1)) {
3473          CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
3474          CoinMemcpyN(savedSolution_+numberColumns_ ,
3475                  numberRows_,rowActivityWork_);
3476          CoinMemcpyN(savedSolution_ ,
3477                  numberColumns_,columnActivityWork_);
3478          // restore extra stuff
3479          int dummy;
3480          matrix_->generalExpanded(this,6,dummy);
3481          // debug
3482          int returnCode = internalFactorize(1);
3483          while (returnCode) {
3484            // ouch
3485            // switch off dense
3486            int saveDense = factorization_->denseThreshold();
3487            factorization_->setDenseThreshold(0);
3488            // Go to safe
3489            factorization_->pivotTolerance(0.99);
3490            // make sure will do safe factorization
3491            pivotVariable_[0]=-1;
3492            returnCode=internalFactorize(2);
3493            factorization_->setDenseThreshold(saveDense);
3494          }
3495        }
3496        type = 2; // so will restore weights
3497        // get primal and dual solutions
3498        gutsOfSolution(givenDuals,NULL);
3499      } 
3500    }
3501#endif
3502  }
3503  // Up tolerance if looks a bit odd
3504  if (numberIterations_>CoinMax(1000,numberRows_>>4)&&(specialOptions_&64)!=0) {
3505    if (sumPrimalInfeasibilities_&&sumPrimalInfeasibilities_<1.0e5) {
3506      int backIteration = progress_->lastIterationNumber(CLP_PROGRESS-1);
3507      if (backIteration>0&&numberIterations_-backIteration<9*CLP_PROGRESS) {
3508        if (factorization_->pivotTolerance()<0.9) {
3509          // up tolerance
3510          factorization_->pivotTolerance(CoinMin(factorization_->pivotTolerance()*1.05+0.02,0.91));
3511          //printf("tol now %g\n",factorization_->pivotTolerance());
3512          progress_->clearIterationNumbers();
3513        }
3514      }
3515    }
3516  }
3517  // Check if looping
3518  int loop;
3519  if (!givenDuals&&type!=2) 
3520    loop = progress_->looping();
3521  else
3522    loop=-1;
3523  int situationChanged=0;
3524  if (loop>=0) {
3525    problemStatus_ = loop; //exit if in loop
3526    if (!problemStatus_) {
3527      // declaring victory
3528      numberPrimalInfeasibilities_ = 0;
3529      sumPrimalInfeasibilities_ = 0.0;
3530    } else {
3531      problemStatus_ = 10; // instead - try other algorithm
3532    }
3533    return;
3534  } else if (loop<-1) {
3535    // something may have changed
3536    gutsOfSolution(NULL,NULL);
3537    situationChanged=1;
3538  }
3539  // really for free variables in
3540  if((progressFlag_&2)!=0) {
3541    situationChanged=2;
3542  }
3543  progressFlag_ = 0; //reset progress flag
3544  if (handler_->detail(CLP_SIMPLEX_STATUS,messages_)<100) {
3545    handler_->message(CLP_SIMPLEX_STATUS,messages_)
3546      <<numberIterations_<<objectiveValue();
3547    handler_->printing(sumPrimalInfeasibilities_>0.0)
3548      <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
3549    handler_->printing(sumDualInfeasibilities_>0.0)
3550      <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
3551    handler_->printing(numberDualInfeasibilitiesWithoutFree_
3552                       <numberDualInfeasibilities_)
3553                         <<numberDualInfeasibilitiesWithoutFree_;
3554    handler_->message()<<CoinMessageEol;
3555  }
3556  realDualInfeasibilities=sumDualInfeasibilities_;
3557  double saveTolerance =dualTolerance_;
3558  // If we need to carry on cleaning variables
3559  if (!numberPrimalInfeasibilities_&&(specialOptions_&1024)!=0&&CLEAN_FIXED) {
3560    for (int iRow=0;iRow<numberRows_;iRow++) {
3561      int iPivot=pivotVariable_[iRow];
3562      if (!flagged(iPivot)&&pivoted(iPivot)) {
3563        // carry on
3564        numberPrimalInfeasibilities_=-1;
3565        sumOfRelaxedPrimalInfeasibilities_ = 1.0;
3566        sumPrimalInfeasibilities_ = 1.0;
3567        break;
3568      }
3569    }
3570  }
3571  /* If we are primal feasible and any dual infeasibilities are on
3572     free variables then it is better to go to primal */
3573  if (!numberPrimalInfeasibilities_&&!numberDualInfeasibilitiesWithoutFree_&&
3574      numberDualInfeasibilities_)
3575    problemStatus_=10;
3576  // dual bound coming in
3577  //double saveDualBound = dualBound_;
3578  while (problemStatus_<=-3) {
3579    int cleanDuals=0;
3580    if (situationChanged!=0)
3581      cleanDuals=1;
3582    int numberChangedBounds=0;
3583    int doOriginalTolerance=0;
3584    if ( lastCleaned==numberIterations_)
3585      doOriginalTolerance=1;
3586    // check optimal
3587    // give code benefit of doubt
3588    if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
3589        sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
3590      // say optimal (with these bounds etc)
3591      numberDualInfeasibilities_ = 0;
3592      sumDualInfeasibilities_ = 0.0;
3593      numberPrimalInfeasibilities_ = 0;
3594      sumPrimalInfeasibilities_ = 0.0;
3595    }
3596    //if (dualFeasible()||problemStatus_==-4||(primalFeasible()&&!numberDualInfeasibilitiesWithoutFree_)) {
3597    if (dualFeasible()||problemStatus_==-4) {
3598      progress_->modifyObjective(objectiveValue_
3599                               -sumDualInfeasibilities_*dualBound_);
3600      if (primalFeasible()&&!givenDuals) {
3601        normalType=false;
3602        // may be optimal - or may be bounds are wrong
3603        handler_->message(CLP_DUAL_BOUNDS,messages_)
3604          <<dualBound_
3605          <<CoinMessageEol;
3606        // save solution in case unbounded
3607        CoinMemcpyN(columnActivityWork_,numberColumns_,
3608                          columnArray_[0]->denseVector());
3609        CoinMemcpyN(rowActivityWork_,numberRows_,
3610                          rowArray_[2]->denseVector());
3611        numberChangedBounds=changeBounds(false,rowArray_[3],changeCost);
3612        if (numberChangedBounds<=0&&!numberDualInfeasibilities_) {
3613          //looks optimal - do we need to reset tolerance
3614          if (perturbation_==101) {
3615            perturbation_=102; // stop any perturbations
3616            cleanDuals=1;
3617            // make sure fake bounds are back
3618            //computeObjectiveValue();
3619            changeBounds(true,NULL,changeCost);
3620            //computeObjectiveValue();
3621            createRim(4);
3622            // make sure duals are current
3623            computeDuals(givenDuals);
3624            checkDualSolution();
3625            //if (numberDualInfeasibilities_)
3626              numberChanged_=1; // force something to happen
3627            //else
3628            //computeObjectiveValue();
3629          }
3630          if (lastCleaned<numberIterations_&&numberTimesOptimal_<4&&
3631              (numberChanged_||(specialOptions_&4096)==0)) {
3632            doOriginalTolerance=2;
3633            numberTimesOptimal_++;
3634            changeMade_++; // say something changed
3635            if (numberTimesOptimal_==1) {
3636              dualTolerance_ = dblParam_[ClpDualTolerance];
3637              // better to have small tolerance even if slower
3638              factorization_->zeroTolerance(1.0e-15);
3639            } else {
3640              dualTolerance_ = dblParam_[ClpDualTolerance];
3641              dualTolerance_ *= pow(2.0,numberTimesOptimal_-1);
3642            }
3643            cleanDuals=2; // If nothing changed optimal else primal
3644          } else {
3645            problemStatus_=0; // optimal
3646            if (lastCleaned<numberIterations_&&numberChanged_) {
3647              handler_->message(CLP_SIMPLEX_GIVINGUP,messages_)
3648                <<CoinMessageEol;
3649            }
3650          }
3651        } else {
3652          cleanDuals=1;
3653          if (doOriginalTolerance==1) {
3654            // check unbounded
3655            // find a variable with bad dj
3656            int iSequence;
3657            int iChosen=-1;
3658            double largest = 100.0*primalTolerance_;
3659            for (iSequence=0;iSequence<numberRows_+numberColumns_;
3660                 iSequence++) {
3661              double djValue = dj_[iSequence];
3662              double originalLo = originalLower(iSequence);
3663              double originalUp = originalUpper(iSequence);
3664              if (fabs(djValue)>fabs(largest)) {
3665                if (getStatus(iSequence)!=basic) {
3666                  if (djValue>0&&originalLo<-1.0e20) {
3667                    if (djValue>fabs(largest)) {
3668                      largest=djValue;
3669                      iChosen=iSequence;
3670                    }
3671                  } else if (djValue<0&&originalUp>1.0e20) {
3672                    if (-djValue>fabs(largest)) {
3673                      largest=djValue;
3674                      iChosen=iSequence;
3675                    }
3676                  } 
3677                }
3678              }
3679            }
3680            if (iChosen>=0) {
3681              int iSave=sequenceIn_;
3682              sequenceIn_=iChosen;
3683              unpack(rowArray_[1]);
3684              sequenceIn_ = iSave;
3685              // if dual infeasibilities then must be free vector so add in dual
3686              if (numberDualInfeasibilities_) {
3687                if (fabs(changeCost)>1.0e-5)
3688                  printf("Odd free/unbounded combo\n");
3689                changeCost += cost_[iChosen];
3690              }
3691              problemStatus_ = checkUnbounded(rowArray_[1],rowArray_[0],
3692                                              changeCost);
3693              rowArray_[1]->clear();
3694            } else {
3695              problemStatus_=-3;
3696            }
3697            if (problemStatus_==2&&perturbation_==101) {
3698              perturbation_=102; // stop any perturbations
3699              cleanDuals=1;
3700              createRim(4);
3701              problemStatus_=-1;
3702            }
3703            if (problemStatus_==2) {
3704              // it is unbounded - restore solution
3705              // but first add in changes to non-basic
3706              int iColumn;
3707              double * original = columnArray_[0]->denseVector();
3708              for (iColumn=0;iColumn<numberColumns_;iColumn++) {
3709                if(getColumnStatus(iColumn)!= basic)
3710                  ray_[iColumn] += 
3711                    columnActivityWork_[iColumn]-original[iColumn];
3712                columnActivityWork_[iColumn] = original[iColumn];
3713              }
3714              CoinMemcpyN(rowArray_[2]->denseVector(),numberRows_,
3715                                rowActivityWork_);
3716            }
3717          } else {
3718            doOriginalTolerance=2;
3719            rowArray_[0]->clear();
3720          }
3721        }
3722        CoinZeroN(columnArray_[0]->denseVector(),numberColumns_);
3723        CoinZeroN(rowArray_[2]->denseVector(),numberRows_);
3724      } 
3725      if (problemStatus_==-4||problemStatus_==-5) {
3726        // may be infeasible - or may be bounds are wrong
3727        numberChangedBounds=changeBounds(false,NULL,changeCost);
3728        /* Should this be here as makes no difference to being feasible.
3729           But seems to make a difference to run times. */
3730        if (perturbation_==101&&0) {
3731          perturbation_=102; // stop any perturbations
3732          cleanDuals=1;
3733          numberChangedBounds=1;
3734          // make sure fake bounds are back
3735          changeBounds(true,NULL,changeCost);
3736          createRim(4);
3737        }
3738        if (numberChangedBounds<=0||dualBound_>1.0e20||
3739            (largestPrimalError_>1.0&&dualBound_>1.0e17)) {
3740          problemStatus_=1; // infeasible
3741          if (perturbation_==101) {
3742            perturbation_=102; // stop any perturbations
3743            //cleanDuals=1;
3744            //numberChangedBounds=1;
3745            //createRim(4);
3746          }
3747        } else {
3748          normalType=false;
3749          problemStatus_=-1; //iterate
3750          cleanDuals=1;
3751          if (numberChangedBounds<=0)
3752            doOriginalTolerance=2;
3753          // and delete ray which has been created
3754          delete [] ray_;
3755          ray_ = NULL;
3756        }
3757
3758      }
3759    } else {
3760      cleanDuals=1;
3761    }
3762    if (problemStatus_<0) {
3763      if (doOriginalTolerance==2) {
3764        // put back original tolerance
3765        lastCleaned=numberIterations_;
3766        numberChanged_ =0; // Number of variables with changed costs
3767        handler_->message(CLP_DUAL_ORIGINAL,messages_)
3768          <<CoinMessageEol;
3769        perturbation_=102; // stop any perturbations
3770#if 0
3771        double * xcost = new double[numberRows_+numberColumns_];
3772        double * xlower = new double[numberRows_+numberColumns_];
3773        double * xupper = new double[numberRows_+numberColumns_];
3774        double * xdj = new double[numberRows_+numberColumns_];
3775        double * xsolution = new double[numberRows_+numberColumns_];
3776        memcpy(xcost,cost_,(numberRows_+numberColumns_)*sizeof(double));
3777        memcpy(xlower,lower_,(numberRows_+numberColumns_)*sizeof(double));
3778        memcpy(xupper,upper_,(numberRows_+numberColumns_)*sizeof(double));
3779        memcpy(xdj,dj_,(numberRows_+numberColumns_)*sizeof(double));
3780        memcpy(xsolution,solution_,(numberRows_+numberColumns_)*sizeof(double));
3781#endif
3782        createRim(4);
3783        // make sure duals are current
3784        computeDuals(givenDuals);
3785        checkDualSolution();
3786#if 0
3787        int i;
3788        for (i=0;i<numberRows_+numberColumns_;i++) {
3789          if (cost_[i]!=xcost[i])
3790            printf("** %d old cost %g new %g sol %g\n",
3791                   i,xcost[i],cost_[i],solution_[i]);
3792          if (lower_[i]!=xlower[i])
3793            printf("** %d old lower %g new %g sol %g\n",
3794                   i,xlower[i],lower_[i],solution_[i]);
3795          if (upper_[i]!=xupper[i])
3796            printf("** %d old upper %g new %g sol %g\n",
3797                   i,xupper[i],upper_[i],solution_[i]);
3798          if (dj_[i]!=xdj[i])
3799            printf("** %d old dj %g new %g sol %g\n",
3800                   i,xdj[i],dj_[i],solution_[i]);
3801          if (solution_[i]!=xsolution[i])
3802            printf("** %d old solution %g new %g sol %g\n",
3803                   i,xsolution[i],solution_[i],solution_[i]);
3804        }
3805        //delete [] xcost;
3806        //delete [] xupper;
3807        //delete [] xlower;
3808        //delete [] xdj;
3809        //delete [] xsolution;
3810#endif
3811        // put back bounds as they were if was optimal
3812        if (doOriginalTolerance==2&&cleanDuals!=2) {
3813          changeMade_++; // say something changed
3814          /* We may have already changed some bounds in this function
3815             so save numberFake_ and add in.
3816
3817             Worst that can happen is that we waste a bit of time  - but it must be finite.
3818          */
3819          int saveNumberFake = numberFake_;
3820          changeBounds(true,NULL,changeCost);
3821          numberFake_ += saveNumberFake;
3822          cleanDuals=2;
3823          //cleanDuals=1;
3824        }
3825#if 0
3826        //int i;
3827        for (i=0;i<numberRows_+numberColumns_;i++) {
3828          if (cost_[i]!=xcost[i])
3829            printf("** %d old cost %g new %g sol %g\n",
3830                   i,xcost[i],cost_[i],solution_[i]);
3831          if (lower_[i]!=xlower[i])
3832            printf("** %d old lower %g new %g sol %g\n",
3833                   i,xlower[i],lower_[i],solution_[i]);
3834          if (upper_[i]!=xupper[i])
3835            printf("** %d old upper %g new %g sol %g\n",
3836                   i,xupper[i],upper_[i],solution_[i]);
3837          if (dj_[i]!=xdj[i])
3838            printf("** %d old dj %g new %g sol %g\n",
3839                   i,xdj[i],dj_[i],solution_[i]);
3840          if (solution_[i]!=xsolution[i])
3841            printf("** %d old solution %g new %g sol %g\n",
3842                   i,xsolution[i],solution_[i],solution_[i]);
3843        }
3844        delete [] xcost;
3845        delete [] xupper;
3846        delete [] xlower;
3847        delete [] xdj;
3848        delete [] xsolution;
3849#endif
3850      }
3851      if (cleanDuals==1||(cleanDuals==2&&!numberDualInfeasibilities_)) {
3852        // make sure dual feasible
3853        // look at all rows and columns
3854        rowArray_[0]->clear();
3855        columnArray_[0]->clear();
3856        double objectiveChange=0.0;
3857#if 0
3858        double * xcost = new double[numberRows_+numberColumns_];
3859        double * xlower = new double[numberRows_+numberColumns_];
3860        double * xupper = new double[numberRows_+numberColumns_];
3861        double * xdj = new double[numberRows_+numberColumns_];
3862        double * xsolution = new double[numberRows_+numberColumns_];
3863        memcpy(xcost,cost_,(numberRows_+numberColumns_)*sizeof(double));
3864        memcpy(xlower,lower_,(numberRows_+numberColumns_)*sizeof(double));
3865        memcpy(xupper,upper_,(numberRows_+numberColumns_)*sizeof(double));
3866        memcpy(xdj,dj_,(numberRows_+numberColumns_)*sizeof(double));
3867        memcpy(xsolution,solution_,(numberRows_+numberColumns_)*sizeof(double));
3868#endif
3869        if (givenDuals)
3870          dualTolerance_=1.0e50;
3871        updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[1],
3872          0.0,objectiveChange,true);
3873        dualTolerance_=saveTolerance;
3874#if 0
3875        int i;
3876        for (i=0;i<numberRows_+numberColumns_;i++) {
3877          if (cost_[i]!=xcost[i])
3878            printf("** %d old cost %g new %g sol %g\n",
3879                   i,xcost[i],cost_[i],solution_[i]);
3880          if (lower_[i]!=xlower[i])
3881            printf("** %d old lower %g new %g sol %g\n",
3882                   i,xlower[i],lower_[i],solution_[i]);
3883          if (upper_[i]!=xupper[i])
3884            printf("** %d old upper %g new %g sol %g\n",
3885                   i,xupper[i],upper_[i],solution_[i]);
3886          if (dj_[i]!=xdj[i])
3887            printf("** %d old dj %g new %g sol %g\n",
3888                   i,xdj[i],dj_[i],solution_[i]);
3889          if (solution_[i]!=xsolution[i])
3890            printf("** %d old solution %g new %g sol %g\n",
3891                   i,xsolution[i],solution_[i],solution_[i]);
3892        }
3893        delete [] xcost;
3894        delete [] xupper;
3895        delete [] xlower;
3896        delete [] xdj;
3897        delete [] xsolution;
3898#endif
3899        // for now - recompute all
3900        gutsOfSolution(NULL,NULL);
3901        if (givenDuals)
3902          dualTolerance_=1.0e50;
3903        updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[1],
3904          0.0,objectiveChange,true);
3905        dualTolerance_=saveTolerance;
3906        //assert(numberDualInfeasibilitiesWithoutFree_==0);
3907
3908        if (numberDualInfeasibilities_||situationChanged==2) 
3909          problemStatus_=-1; // carry on as normal
3910        situationChanged=0;
3911      } else {
3912        // iterate
3913        if (cleanDuals!=2) 
3914          problemStatus_=-1;
3915        else 
3916          problemStatus_ = 10; // try primal
3917      }
3918    }
3919  }
3920  if (type==0||type==1) {
3921    if (!type) {
3922      // create save arrays
3923      delete [] saveStatus_;
3924      delete [] savedSolution_;
3925      saveStatus_ = new unsigned char [numberRows_+numberColumns_];
3926      savedSolution_ = new double [numberRows_+numberColumns_];
3927    }
3928    // save arrays
3929    CoinMemcpyN(status_,numberColumns_+numberRows_,saveStatus_);
3930    CoinMemcpyN(rowActivityWork_,
3931           numberRows_,savedSolution_+numberColumns_);
3932    CoinMemcpyN(columnActivityWork_,numberColumns_,savedSolution_);
3933    // save extra stuff
3934    int dummy;
3935    matrix_->generalExpanded(this,5,dummy);
3936  }
3937
3938  // restore weights (if saved) - also recompute infeasibility list
3939  if (tentativeStatus>-3) 
3940    dualRowPivot_->saveWeights(this,(type <2) ? 2 : 4);
3941  else
3942    dualRowPivot_->saveWeights(this,3);
3943  // unflag all variables (we may want to wait a bit?)
3944  if ((tentativeStatus!=-2&&tentativeStatus!=-1)&&unflagVariables) {
3945    int iRow;
3946    int numberFlagged=0;
3947    for (iRow=0;iRow<numberRows_;iRow++) {
3948      int iPivot=pivotVariable_[iRow];
3949      if (flagged(iPivot)) {
3950        numberFlagged++;
3951        clearFlagged(iPivot);
3952      }
3953    }
3954#if 0
3955    if (numberFlagged) {
3956      printf("unflagging %d variables - status %d ninf %d nopt %d\n",numberFlagged,tentativeStatus,
3957             numberPrimalInfeasibilities_,
3958             numberTimesOptimal_);
3959    }
3960#endif
3961    unflagVariables = numberFlagged>0;
3962    if (numberFlagged&&!numberPivots) {
3963      /* looks like trouble as we have not done any iterations.
3964         Try changing pivot tolerance then give it a few goes and give up */
3965      if (factorization_->pivotTolerance()<0.9) {
3966        factorization_->pivotTolerance(0.99);
3967      } else if (numberTimesOptimal_<10) {
3968        numberTimesOptimal_++;
3969      } else {
3970        unflagVariables=false;
3971        changeMade_=false;
3972        secondaryStatus_ = 1; // and say probably infeasible
3973      }
3974    }
3975  }
3976  // see if cutoff reached
3977  double limit = 0.0;
3978  getDblParam(ClpDualObjectiveLimit, limit);
3979  if(fabs(limit)<1.0e30&&objectiveValue()*optimizationDirection_>
3980           limit&&
3981           !numberAtFakeBound()&&!numberDualInfeasibilities_) {
3982    //printf("lim %g obj %g %g\n",limit,objectiveValue_,objectiveValue());
3983    if (true||perturbation_==101) {
3984      // be careful
3985      if (numberIterations_) {
3986        // I don't want to change calling sequences - so ...
3987        int iSequence;
3988        objectiveValue_ = 0.0;
3989        const double * obj = objective();
3990        if (!columnScale_) {
3991          for (iSequence=0;iSequence<numberColumns_;iSequence++) {
3992            double value = columnActivityWork_[iSequence];
3993            objectiveValue_ += value*obj[iSequence];
3994          }
3995        } else {
3996          for (iSequence=0;iSequence<numberColumns_;iSequence++) {
3997            double scaleFactor = columnScale_[iSequence];
3998            double valueScaled = columnActivityWork_[iSequence];
3999            objectiveValue_ += valueScaled*scaleFactor*obj[iSequence];
4000          }
4001        }
4002        // But remember direction as we are using external objective
4003        objectiveValue_ *= optimizationDirection_;
4004        objectiveValue_ += objective_->nonlinearOffset();
4005        objectiveValue_ /= (objectiveScale_*rhsScale_);
4006        //computeObjectiveValue(); // value without perturbation
4007        if(objectiveValue()*optimizationDirection_>limit) {
4008          problemStatus_=1;
4009          secondaryStatus_ = 1; // and say was on cutoff
4010        }
4011      }
4012    } else {
4013      problemStatus_=1;
4014      secondaryStatus_ = 1; // and say was on cutoff
4015    }
4016  }
4017  // If we are in trouble and in branch and bound give up
4018  if ((specialOptions_&1024)!=0) {
4019    int looksBad=0;
4020    if (largestPrimalError_*largestDualError_>1.0e2) {
4021      looksBad=1;
4022    } else if (largestPrimalError_>1.0e-2
4023        &&objectiveValue_>CoinMin(1.0e15,1.0e3*limit)) {
4024      looksBad=2;
4025    }
4026    if (looksBad) {
4027      if (factorization_->pivotTolerance()<0.9) {
4028        // up tolerance
4029        factorization_->pivotTolerance(CoinMin(factorization_->pivotTolerance()*1.05+0.02,0.91));
4030      } else if (numberIterations_>10000) {
4031        if (handler_->logLevel()>0)
4032          printf("bad dual - saying infeasible %d\n",looksBad);
4033        problemStatus_=1;
4034        secondaryStatus_ = 1; // and say was on cutoff
4035      } else if (largestPrimalError_>1.0e5) {
4036        allSlackBasis();
4037        problemStatus_=10;
4038        //if (handler_->logLevel()>0)
4039          printf("bad dual - going to primal %d %g\n",looksBad,largestPrimalError_);
4040      }
4041    }
4042  }
4043  if (problemStatus_<0&&!changeMade_) {
4044    problemStatus_=4; // unknown
4045  }
4046  lastGoodIteration_ = numberIterations_;
4047  if (problemStatus_<0) {
4048    sumDualInfeasibilities_=realDualInfeasibilities; // back to say be careful
4049    if (sumDualInfeasibilities_)
4050      numberDualInfeasibilities_=1;
4051  }
4052#if 1
4053  double thisObj = progress_->lastObjective(0);
4054  double lastObj = progress_->lastObjective(1);
4055  if (lastObj>thisObj+1.0e-4*CoinMax(fabs(thisObj),fabs(lastObj))+1.0e-4
4056      &&givenDuals==NULL&&firstFree_<0) {
4057    int maxFactor = factorization_->maximumPivots();
4058    if (maxFactor>10) {
4059      if (forceFactorization_<0)
4060        forceFactorization_= maxFactor;
4061      forceFactorization_ = CoinMax(1,(forceFactorization_>>1));
4062      //printf("Reducing factorization frequency\n");
4063    } 
4064  }
4065#endif
4066}
4067/* While updateDualsInDual sees what effect is of flip
4068   this does actual flipping.
4069   If change >0.0 then value in array >0.0 => from lower to upper
4070*/
4071void 
4072ClpSimplexDual::flipBounds(CoinIndexedVector * rowArray,
4073                  CoinIndexedVector * columnArray,
4074                  double change)
4075{
4076  int number;
4077  int * which;
4078 
4079  int iSection;
4080
4081  for (iSection=0;iSection<2;iSection++) {
4082    int i;
4083    double * solution = solutionRegion(iSection);
4084    double * lower = lowerRegion(iSection);
4085    double * upper = upperRegion(iSection);
4086    int addSequence;
4087    if (!iSection) {
4088      number = rowArray->getNumElements();
4089      which = rowArray->getIndices();
4090      addSequence = numberColumns_;
4091    } else {
4092      number = columnArray->getNumElements();
4093      which = columnArray->getIndices();
4094      addSequence = 0;
4095    }
4096   
4097    for (i=0;i<number;i++) {
4098      int iSequence = which[i];
4099      Status status = getStatus(iSequence+addSequence);
4100
4101      switch(status) {
4102
4103      case basic:
4104      case isFree:
4105      case superBasic:
4106      case ClpSimplex::isFixed:
4107        break;
4108      case atUpperBound:
4109        // to lower bound
4110        setStatus(iSequence+addSequence,atLowerBound);
4111        solution[iSequence] = lower[iSequence];
4112        break;
4113      case atLowerBound:
4114        // to upper bound
4115        setStatus(iSequence+addSequence,atUpperBound);
4116        solution[iSequence] = upper[iSequence];
4117        break;
4118      }
4119    }
4120  }
4121  rowArray->setNumElements(0);
4122  columnArray->setNumElements(0);
4123}
4124// Restores bound to original bound
4125void 
4126ClpSimplexDual::originalBound( int iSequence)
4127{
4128  if (getFakeBound(iSequence)!=noFake)
4129    numberFake_--;;
4130  if (auxiliaryModel_) {
4131    // just copy back
4132    lower_[iSequence]=auxiliaryModel_->lowerRegion()[iSequence+numberRows_+numberColumns_];
4133    upper_[iSequence]=auxiliaryModel_->upperRegion()[iSequence+numberRows_+numberColumns_];
4134    return;
4135  }
4136  if (iSequence>=numberColumns_) {
4137    // rows
4138    int iRow = iSequence-numberColumns_;
4139    rowLowerWork_[iRow]=rowLower_[iRow];
4140    rowUpperWork_[iRow]=rowUpper_[iRow];
4141    if (rowScale_) {
4142      if (rowLowerWork_[iRow]>-1.0e50)
4143        rowLowerWork_[iRow] *= rowScale_[iRow]*rhsScale_;
4144      if (rowUpperWork_[iRow]<1.0e50)
4145        rowUpperWork_[iRow] *= rowScale_[iRow]*rhsScale_;
4146    } else if (rhsScale_!=1.0) {
4147      if (rowLowerWork_[iRow]>-1.0e50)
4148        rowLowerWork_[iRow] *= rhsScale_;
4149      if (rowUpperWork_[iRow]<1.0e50)
4150        rowUpperWork_[iRow] *= rhsScale_;
4151    }
4152  } else {
4153    // columns
4154    columnLowerWork_[iSequence]=columnLower_[iSequence];
4155    columnUpperWork_[iSequence]=columnUpper_[iSequence];
4156    if (rowScale_) {
4157      double multiplier = 1.0/columnScale_[iSequence];
4158      if (columnLowerWork_[iSequence]>-1.0e50)
4159        columnLowerWork_[iSequence] *= multiplier*rhsScale_;
4160      if (columnUpperWork_[iSequence]<1.0e50)
4161        columnUpperWork_[iSequence] *= multiplier*rhsScale_;
4162    } else if (rhsScale_!=1.0) {
4163      if (columnLowerWork_[iSequence]>-1.0e50)
4164        columnLowerWork_[iSequence] *= rhsScale_;
4165      if (columnUpperWork_[iSequence]<1.0e50)
4166        columnUpperWork_[iSequence] *= rhsScale_;
4167    }
4168  }
4169  setFakeBound(iSequence,noFake);
4170}
4171/* As changeBounds but just changes new bounds for a single variable.
4172   Returns true if change */
4173bool 
4174ClpSimplexDual::changeBound( int iSequence)
4175{
4176  // old values
4177  double oldLower=lower_[iSequence];
4178  double oldUpper=upper_[iSequence];
4179  double value=solution_[iSequence];
4180  bool modified=false;
4181  originalBound(iSequence);
4182  // original values
4183  double lowerValue=lower_[iSequence];
4184  double upperValue=upper_[iSequence];
4185  // back to altered values
4186  lower_[iSequence] = oldLower;
4187  upper_[iSequence] = oldUpper;
4188  if (getFakeBound(iSequence)!=noFake)
4189    numberFake_--;;
4190  if (value==oldLower) {
4191    if (upperValue > oldLower + dualBound_) {
4192      upper_[iSequence]=oldLower+dualBound_;
4193      setFakeBound(iSequence,upperFake);
4194      modified=true;
4195      numberFake_++;
4196    }
4197  } else if (value==oldUpper) {
4198    if (lowerValue < oldUpper - dualBound_) {
4199      lower_[iSequence]=oldUpper-dualBound_;
4200      setFakeBound(iSequence,lowerFake);
4201      modified=true;
4202      numberFake_++;
4203    }
4204  } else {
4205    assert(value==oldLower||value==oldUpper);
4206  }
4207  return modified;
4208}
4209// Perturbs problem
4210void 
4211ClpSimplexDual::perturb()
4212{
4213  if (perturbation_>100)
4214    return; //perturbed already
4215  if (perturbation_==100)
4216    perturbation_=50; // treat as normal
4217  int savePerturbation = perturbation_;
4218  bool modifyRowCosts=false;
4219  // dual perturbation
4220  double perturbation=1.0e-20;
4221  // maximum fraction of cost to perturb
4222  double maximumFraction = 1.0e-5;
4223  double constantPerturbation = 100.0*dualTolerance_;
4224  int maxLength=0;
4225  int minLength=numberRows_;
4226  double averageCost = 0.0;
4227  // look at element range
4228  double smallestNegative;
4229  double largestNegative;
4230  double smallestPositive;
4231  double largestPositive;
4232  matrix_->rangeOfElements(smallestNegative, largestNegative,
4233                           smallestPositive, largestPositive);
4234  smallestPositive = CoinMin(fabs(smallestNegative),smallestPositive);
4235  largestPositive = CoinMax(fabs(largestNegative),largestPositive);
4236  double elementRatio = largestPositive/smallestPositive;
4237  int numberNonZero=0;
4238  if (!numberIterations_&&perturbation_==50) {
4239    // See if we need to perturb
4240    double * sort = new double[numberColumns_];
4241    // Use objective BEFORE scaling
4242    const double * obj = objective();
4243    int i;
4244    for (i=0;i<numberColumns_;i++) {
4245      double value = fabs(obj[i]);
4246      sort[i]=value;
4247      averageCost += value;
4248      if (value)
4249        numberNonZero++;
4250    }
4251    if (numberNonZero)
4252      averageCost /= (double) numberNonZero;
4253    else
4254      averageCost = 1.0;
4255    std::sort(sort,sort+numberColumns_);
4256    int number=1;
4257    double last = sort[0];
4258    for (i=1;i<numberColumns_;i++) {
4259      if (last!=sort[i])
4260        number++;
4261      last=sort[i];
4262    }
4263    delete [] sort;
4264#if 0
4265    printf("nnz %d percent %d",number,(number*100)/numberColumns_);
4266    if (number*4>numberColumns_)
4267      printf(" - Would not perturb\n");
4268    else
4269      printf(" - Would perturb\n");
4270    //exit(0);
4271#endif
4272    //printf("ratio number diff costs %g, element ratio %g\n",((double)number)/((double) numberColumns_),
4273    //                                                                elementRatio);
4274    //number=0;
4275    if (number*4>numberColumns_||elementRatio>1.0e12) {
4276      perturbation_=100;
4277      return; // good enough
4278    }
4279  }
4280  int iColumn;
4281  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4282    if (columnLowerWork_[iColumn]<columnUpperWork_[iColumn]) {
4283      int length = matrix_->getVectorLength(iColumn);
4284      if (length>2) {
4285        maxLength = CoinMax(maxLength,length);
4286        minLength = CoinMin(minLength,length);
4287      }
4288    }
4289  }
4290  // If > 70 then do rows
4291  if (perturbation_>=70) {
4292    modifyRowCosts=true;
4293    perturbation_ -= 20;
4294    printf("Row costs modified, ");
4295  }
4296  bool uniformChange=false;
4297  if (perturbation_>50) {
4298    // Experiment
4299    // maximumFraction could be 1.0e-10 to 1.0
4300    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};
4301    maximumFraction = m[CoinMin(perturbation_-51,10)];
4302  }
4303  int iRow;
4304  double smallestNonZero=1.0e100;
4305  numberNonZero=0;
4306  if (perturbation_>=50) {
4307    perturbation = 1.0e-8;
4308    bool allSame=true;
4309    double lastValue=0.0;
4310    for (iRow=0;iRow<numberRows_;iRow++) {
4311      double lo = rowLowerWork_[iRow];
4312      double up = rowUpperWork_[iRow];
4313      if (lo<up) {
4314        double value = fabs(rowObjectiveWork_[iRow]);
4315        perturbation = CoinMax(perturbation,value);
4316        if (value) {
4317          modifyRowCosts=true;
4318          smallestNonZero = CoinMin(smallestNonZero,value);
4319        }
4320      } 
4321      if (lo&&lo>-1.0e10) {
4322        numberNonZero++;
4323        lo=fabs(lo);
4324        if (!lastValue) 
4325          lastValue=lo;
4326        else if (fabs(lo-lastValue)>1.0e-7)
4327          allSame=false;
4328      }
4329      if (up&&up<1.0e10) {
4330        numberNonZero++;
4331        up=fabs(up);
4332        if (!lastValue) 
4333          lastValue=up;
4334        else if (fabs(up-lastValue)>1.0e-7)
4335          allSame=false;
4336      }
4337    }
4338    double lastValue2=0.0;
4339    for (iColumn=0;iColumn<numberColumns_;iColumn++) { 
4340      double lo = columnLowerWork_[iColumn];
4341      double up = columnUpperWork_[iColumn];
4342      if (lo<up) {
4343        double value = 
4344          fabs(objectiveWork_[iColumn]);
4345        perturbation = CoinMax(perturbation,value);
4346        if (value) {
4347          smallestNonZero = CoinMin(smallestNonZero,value);
4348        }
4349      }
4350      if (lo&&lo>-1.0e10) {
4351        //numberNonZero++;
4352        lo=fabs(lo);
4353        if (!lastValue2) 
4354          lastValue2=lo;
4355        else if (fabs(lo-lastValue2)>1.0e-7)
4356          allSame=false;
4357      }
4358      if (up&&up<1.0e10) {
4359        //numberNonZero++;
4360        up=fabs(up);
4361        if (!lastValue2) 
4362          lastValue2=up;
4363        else if (fabs(up-lastValue2)>1.0e-7)
4364          allSame=false;
4365      }
4366    }
4367    if (allSame) {
4368      // Check elements
4369      double smallestNegative;
4370      double largestNegative;
4371      double smallestPositive;
4372      double largestPositive;
4373      matrix_->rangeOfElements(smallestNegative,largestNegative,
4374                               smallestPositive,largestPositive);
4375      if (smallestNegative==largestNegative&&
4376          smallestPositive==largestPositive) {
4377        // Really hit perturbation
4378        double adjust = CoinMin(100.0*maximumFraction,1.0e-3*CoinMax(lastValue,lastValue2));
4379        maximumFraction = CoinMax(adjust,maximumFraction);
4380      }
4381    }
4382    perturbation = CoinMin(perturbation,smallestNonZero/maximumFraction);
4383  } else {
4384    // user is in charge
4385    maximumFraction = 1.0e-1;
4386    // but some experiments
4387    if (perturbation_<=-900) {
4388      modifyRowCosts=true;
4389      perturbation_ += 1000;
4390      printf("Row costs modified, ");
4391    }
4392    if (perturbation_<=-10) {
4393      perturbation_ += 10; 
4394      maximumFraction = 1.0;
4395      if ((-perturbation_)%100>=10) {
4396        uniformChange=true;
4397        perturbation_+=20;
4398      }
4399      while (perturbation_<-10) {
4400        perturbation_ += 100;
4401        maximumFraction *= 1.0e-1;
4402      }
4403    }
4404    perturbation = pow(10.0,perturbation_);
4405  }
4406  double largestZero=0.0;
4407  double largest=0.0;
4408  double largestPerCent=0.0;
4409  // modify costs
4410  bool printOut=(handler_->logLevel()==63);
4411  printOut=false;
4412  if (modifyRowCosts) {
4413    for (iRow=0;iRow<numberRows_;iRow++) {
4414      if (rowLowerWork_[iRow]<rowUpperWork_[iRow]) {
4415        double value = perturbation;
4416        double currentValue = rowObjectiveWork_[iRow];
4417        value = CoinMin(value,maximumFraction*(fabs(currentValue)+1.0e-1*perturbation+1.0e-3));
4418        if (rowLowerWork_[iRow]>-largeValue_) {
4419          if (fabs(rowLowerWork_[iRow])<fabs(rowUpperWork_[iRow])) 
4420            value *= CoinDrand48();
4421          else
4422            value *= -CoinDrand48();
4423        } else if (rowUpperWork_[iRow]<largeValue_) {
4424          value *= -CoinDrand48();
4425        } else {
4426          value=0.0;
4427        }
4428        if (currentValue) {
4429          largest = CoinMax(largest,fabs(value));
4430          if (fabs(value)>fabs(currentValue)*largestPerCent)
4431            largestPerCent=fabs(value/currentValue);
4432        } else {
4433          largestZero=CoinMax(largestZero,fabs(value));
4434        }
4435        if (printOut)
4436          printf("row %d cost %g change %g\n",iRow,rowObjectiveWork_[iRow],value);
4437        rowObjectiveWork_[iRow] += value;
4438      }
4439    }
4440  }
4441  // 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};
4442  // 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};
4443  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};
4444  //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};
4445  double extraWeight=10.0;
4446  // Scale back if wanted
4447  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};
4448  if (constantPerturbation<99.0*dualTolerance_) {
4449    perturbation *= 0.1;
4450    extraWeight=0.5;
4451    memcpy(weight,weight2,sizeof(weight2));
4452  }
4453  // adjust weights if all columns long
4454  double factor=1.0;
4455  if (maxLength) {
4456    factor = 3.0/(double) minLength;
4457  }
4458  // Make variables with more elements more expensive
4459  const double m1 = 0.5;
4460  double smallestAllowed = CoinMin(1.0e-2*dualTolerance_,maximumFraction);
4461  //double largestAllowed = CoinMax(1.0e3*dualTolerance_,maximumFraction*10.0*averageCost);
4462  double largestAllowed = CoinMax(1.0e3*dualTolerance_,maximumFraction*averageCost);
4463  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4464    if (columnLowerWork_[iColumn]<columnUpperWork_[iColumn]&&getStatus(iColumn)!=basic) {
4465      double value = perturbation;
4466      double currentValue = objectiveWork_[iColumn];
4467      value = CoinMin(value,constantPerturbation+maximumFraction*(fabs(currentValue)+1.0e-1*perturbation+1.0e-8));
4468      //value = CoinMin(value,constantPerturbation;+maximumFraction*fabs(currentValue));
4469      double value2 = constantPerturbation+1.0e-1*smallestNonZero;
4470      if (uniformChange) {
4471        value = maximumFraction;
4472        value2=maximumFraction;
4473      }
4474      if (columnLowerWork_[iColumn]>-largeValue_) {
4475        if (fabs(columnLowerWork_[iColumn])<
4476            fabs(columnUpperWork_[iColumn])) {
4477          value *= (1.0-m1 +m1*CoinDrand48());
4478          value2 *= (1.0-m1 +m1*CoinDrand48());
4479        } else {
4480          value *= -(1.0-m1+m1*CoinDrand48());
4481          value2 *= -(1.0-m1+m1*CoinDrand48());
4482        }
4483      } else if (columnUpperWork_[iColumn]<largeValue_) {
4484        value *= -(1.0-m1+m1*CoinDrand48());
4485        value2 *= -(1.0-m1+m1*CoinDrand48());
4486      } else {
4487        value=0.0;
4488      }
4489      if (value) {
4490        int length = matrix_->getVectorLength(iColumn);
4491        if (length>3) {
4492          length = (int) ((double) length * factor);
4493          length = CoinMax(3,length);
4494        }
4495        double multiplier;
4496#if 1
4497        if (length<10)
4498          multiplier=weight[length];
4499        else
4500          multiplier = weight[10];
4501#else
4502        if (length<10)
4503          multiplier=weight[length];
4504        else
4505          multiplier = weight[10]+extraWeight*(length-10);
4506        multiplier *= 0.5;
4507#endif
4508        value *= multiplier;
4509        value = CoinMin(value,value2);
4510        if (savePerturbation<50||savePerturbation>60) {
4511          if (fabs(value)<=dualTolerance_)
4512            value=0.0;
4513        } else if (value) {
4514          // get in range
4515          if (fabs(value)<=smallestAllowed) {
4516            value *= 10.0;
4517            while (fabs(value)<=smallestAllowed) 
4518              value *= 10.0;
4519          } else if (fabs(value)>largestAllowed) {
4520            value *= 0.1;
4521            while (fabs(value)>largestAllowed) 
4522              value *= 0.1;
4523          }
4524        }
4525        if (currentValue) {
4526          largest = CoinMax(largest,fabs(value));
4527          if (fabs(value)>fabs(currentValue)*largestPerCent)
4528            largestPerCent=fabs(value/currentValue);
4529        } else {
4530          largestZero=CoinMax(largestZero,fabs(value));
4531        }
4532        if (printOut)
4533          printf("col %d cost %g change %g\n",iColumn,objectiveWork_[iColumn],value);
4534        objectiveWork_[iColumn] += value;
4535      }
4536    }
4537  }
4538  handler_->message(CLP_SIMPLEX_PERTURB,messages_)
4539    <<100.0*maximumFraction<<perturbation<<largest<<100.0*largestPerCent<<largestZero
4540    <<CoinMessageEol;
4541  // and zero changes
4542  //int nTotal = numberRows_+numberColumns_;
4543  //CoinZeroN(cost_+nTotal,nTotal);
4544  // say perturbed
4545  perturbation_=101;
4546
4547}
4548/* For strong branching.  On input lower and upper are new bounds
4549   while on output they are change in objective function values
4550   (>1.0e50 infeasible).
4551   Return code is 0 if nothing interesting, -1 if infeasible both
4552   ways and +1 if infeasible one way (check values to see which one(s))
4553   Returns -2 if bad factorization
4554*/
4555int ClpSimplexDual::strongBranching(int numberVariables,const int * variables,
4556                                    double * newLower, double * newUpper,
4557                                    double ** outputSolution,
4558                                    int * outputStatus, int * outputIterations,
4559                                    bool stopOnFirstInfeasible,
4560                                    bool alwaysFinish,
4561                                    int startFinishOptions)
4562{
4563  int i;
4564  int returnCode=0;
4565  double saveObjectiveValue = objectiveValue_;
4566  algorithm_ = -1;
4567
4568  //scaling(false);
4569
4570  // put in standard form (and make row copy)
4571  // create modifiable copies of model rim and do optional scaling
4572  createRim(7+8+16+32,true,startFinishOptions);
4573
4574  // change newLower and newUpper if scaled
4575
4576  // Do initial factorization
4577  // and set certain stuff
4578  // We can either set increasing rows so ...IsBasic gives pivot row
4579  // or we can just increment iBasic one by one
4580  // for now let ...iBasic give pivot row
4581  int useFactorization=false;
4582  if ((startFinishOptions&2)!=0&&(whatsChanged_&(2+512))==2+512)
4583    useFactorization=true; // Keep factorization if possible
4584  // switch off factorization if bad
4585  if (pivotVariable_[0]<0)
4586    useFactorization=false;
4587  if (!useFactorization||factorization_->numberRows()!=numberRows_) {
4588    useFactorization = false;
4589    factorization_->increasingRows(2);
4590    // row activities have negative sign
4591    factorization_->slackValue(-1.0);
4592    factorization_->zeroTolerance(1.0e-13);
4593
4594    int factorizationStatus = internalFactorize(0);
4595    if (factorizationStatus<0) {
4596      // some error
4597      // we should either debug or ignore
4598#ifndef NDEBUG
4599      printf("***** ClpDual strong branching factorization error - debug\n");
4600#endif
4601      return -2;
4602    } else if (factorizationStatus&&factorizationStatus<=numberRows_) {
4603      handler_->message(CLP_SINGULARITIES,messages_)
4604        <<factorizationStatus
4605        <<CoinMessageEol;
4606    }
4607  }
4608  // save stuff
4609  ClpFactorization saveFactorization(*factorization_);
4610  // save basis and solution
4611  double * saveSolution = new double[numberRows_+numberColumns_];
4612  CoinMemcpyN(solution_,
4613         numberRows_+numberColumns_,saveSolution);
4614  unsigned char * saveStatus = 
4615    new unsigned char [numberRows_+numberColumns_];
4616  CoinMemcpyN(status_,numberColumns_+numberRows_,saveStatus);
4617  // save bounds as createRim makes clean copies
4618  double * saveLower = new double[numberRows_+numberColumns_];
4619  CoinMemcpyN(lower_,
4620         numberRows_+numberColumns_,saveLower);
4621  double * saveUpper = new double[numberRows_+numberColumns_];
4622  CoinMemcpyN(upper_,
4623         numberRows_+numberColumns_,saveUpper);
4624  double * saveObjective = new double[numberRows_+numberColumns_];
4625  CoinMemcpyN(cost_,
4626         numberRows_+numberColumns_,saveObjective);
4627  int * savePivot = new int [numberRows_];
4628  CoinMemcpyN(pivotVariable_, numberRows_,savePivot);
4629  // need to save/restore weights.
4630
4631  int iSolution = 0;
4632  for (i=0;i<numberVariables;i++) {
4633    int iColumn = variables[i];
4634    double objectiveChange;
4635    double saveBound;
4636   
4637    // try down
4638   
4639    saveBound = columnUpper_[iColumn];
4640    // external view - in case really getting optimal
4641    columnUpper_[iColumn] = newUpper[i];
4642    if (scalingFlag_<=0) 
4643      upper_[iColumn] = newUpper[i]*rhsScale_;
4644    else 
4645      upper_[iColumn] = (newUpper[i]/columnScale_[iColumn])*rhsScale_; // scale
4646    // Start of fast iterations
4647    int status = fastDual(alwaysFinish);
4648    CoinAssert (problemStatus_||objectiveValue_<1.0e50);
4649    // make sure plausible
4650    double obj = CoinMax(objectiveValue_,saveObjectiveValue);
4651    if (status&&problemStatus_!=3) {
4652      // not finished - might be optimal
4653      checkPrimalSolution(rowActivityWork_,columnActivityWork_);
4654      double limit = 0.0;
4655      getDblParam(ClpDualObjectiveLimit, limit);
4656      if (!numberPrimalInfeasibilities_&&obj<limit) { 
4657        problemStatus_=0;
4658      } 
4659      status=problemStatus_;
4660    }
4661    if (status||(problemStatus_==0&&!isDualObjectiveLimitReached())) {
4662      objectiveChange = obj-saveObjectiveValue;
4663    } else {
4664      objectiveChange = 1.0e100;
4665      status=1;
4666    }
4667    if (problemStatus_==3)
4668      status=2;
4669
4670    if (scalingFlag_<=0) {
4671      CoinMemcpyN(solution_,numberColumns_,outputSolution[iSolution]);
4672    } else {
4673      int j;
4674      double * sol = outputSolution[iSolution];
4675      for (j=0;j<numberColumns_;j++) 
4676        sol[j] = solution_[j]*columnScale_[j];
4677    }
4678    outputStatus[iSolution]=status;
4679    outputIterations[iSolution]=numberIterations_;
4680    iSolution++;
4681    // restore
4682    CoinMemcpyN(saveSolution,
4683            numberRows_+numberColumns_,solution_);
4684    CoinMemcpyN(saveStatus,numberColumns_+numberRows_,status_);
4685    CoinMemcpyN(saveLower,
4686            numberRows_+numberColumns_,lower_);
4687    CoinMemcpyN(saveUpper,
4688            numberRows_+numberColumns_,upper_);
4689    CoinMemcpyN(saveObjective,
4690            numberRows_+numberColumns_,cost_);
4691    columnUpper_[iColumn] = saveBound;
4692    CoinMemcpyN(savePivot, numberRows_,pivotVariable_);
4693    delete factorization_;
4694    factorization_ = new ClpFactorization(saveFactorization);
4695
4696    newUpper[i]=objectiveChange;
4697#ifdef CLP_DEBUG
4698    printf("down on %d costs %g\n",iColumn,objectiveChange);
4699#endif
4700         
4701    // try up
4702   
4703    saveBound = columnLower_[iColumn];
4704    // external view - in case really getting optimal
4705    columnLower_[iColumn] = newLower[i];
4706    if (scalingFlag_<=0) 
4707      lower_[iColumn] = newLower[i]*rhsScale_;
4708    else 
4709      lower_[iColumn] = (newLower[i]/columnScale_[iColumn])*rhsScale_; // scale
4710    // Start of fast iterations
4711    status = fastDual(alwaysFinish);
4712    // make sure plausible
4713    obj = CoinMax(objectiveValue_,saveObjectiveValue);
4714    if (status&&problemStatus_!=3) {
4715      // not finished - might be optimal
4716      checkPrimalSolution(rowActivityWork_,columnActivityWork_);
4717      double limit = 0.0;
4718      getDblParam(ClpDualObjectiveLimit, limit);
4719      if (!numberPrimalInfeasibilities_&&obj< limit) { 
4720        problemStatus_=0;
4721      } 
4722      status=problemStatus_;
4723    }
4724    if (status||(problemStatus_==0&&!isDualObjectiveLimitReached())) {
4725      objectiveChange = obj-saveObjectiveValue;
4726    } else {
4727      objectiveChange = 1.0e100;
4728      status=1;
4729    }
4730    if (problemStatus_==3)
4731      status=2;
4732    if (scalingFlag_<=0) {
4733      CoinMemcpyN(solution_,numberColumns_,outputSolution[iSolution]);
4734    } else {
4735      int j;
4736      double * sol = outputSolution[iSolution];
4737      for (j=0;j<numberColumns_;j++) 
4738        sol[j] = solution_[j]*columnScale_[j];
4739    }
4740    outputStatus[iSolution]=status;
4741    outputIterations[iSolution]=numberIterations_;
4742    iSolution++;
4743
4744    // restore
4745    CoinMemcpyN(saveSolution,
4746            numberRows_+numberColumns_,solution_);
4747    CoinMemcpyN(saveStatus,numberColumns_+numberRows_,status_);
4748    CoinMemcpyN(saveLower,
4749            numberRows_+numberColumns_,lower_);
4750    CoinMemcpyN(saveUpper,
4751            numberRows_+numberColumns_,upper_);
4752    CoinMemcpyN(saveObjective,
4753            numberRows_+numberColumns_,cost_);
4754    columnLower_[iColumn] = saveBound;
4755    CoinMemcpyN(savePivot, numberRows_,pivotVariable_);
4756    delete factorization_;
4757    factorization_ = new ClpFactorization(saveFactorization);
4758
4759    newLower[i]=objectiveChange;
4760#ifdef CLP_DEBUG
4761    printf("up on %d costs %g\n",iColumn,objectiveChange);
4762#endif
4763         
4764    /* Possibilities are:
4765       Both sides feasible - store
4766       Neither side feasible - set objective high and exit if desired
4767       One side feasible - change bounds and resolve
4768    */
4769    if (newUpper[i]<1.0e100) {
4770      if(newLower[i]<1.0e100) {
4771        // feasible - no action
4772      } else {
4773        // up feasible, down infeasible
4774        returnCode=1;
4775        if (stopOnFirstInfeasible)
4776          break;
4777      }
4778    } else {
4779      if(newLower[i]<1.0e100) {
4780        // down feasible, up infeasible
4781        returnCode=1;
4782        if (stopOnFirstInfeasible)
4783          break;
4784      } else {
4785        // neither side feasible
4786        returnCode=-1;
4787        break;
4788      }
4789    }
4790  }
4791  delete [] saveSolution;
4792  delete [] saveLower;
4793  delete [] saveUpper;
4794  delete [] saveObjective;
4795  delete [] saveStatus;
4796  delete [] savePivot;
4797  if ((startFinishOptions&1)==0) {
4798    deleteRim(1);
4799    whatsChanged_=0;
4800  } else {
4801    // Original factorization will have been put back by last loop
4802    //delete factorization_;
4803    //factorization_ = new ClpFactorization(saveFactorization);
4804    deleteRim(0);
4805    // mark all as current
4806    whatsChanged_ = 0xffff;
4807  }
4808  objectiveValue_ = saveObjectiveValue;
4809  return returnCode;
4810}
4811// treat no pivot as finished (unless interesting)
4812int ClpSimplexDual::fastDual(bool alwaysFinish)
4813{
4814  algorithm_ = -1;
4815  secondaryStatus_=0;
4816  // Say in fast dual
4817  specialOptions_ |= 16384;
4818  //handler_->setLogLevel(63);
4819  // save data
4820  ClpDataSave data = saveData();
4821  dualTolerance_=dblParam_[ClpDualTolerance];
4822  primalTolerance_=dblParam_[ClpPrimalTolerance];
4823
4824  // save dual bound
4825  double saveDualBound = dualBound_;
4826
4827  double objectiveChange;
4828  // for dual we will change bounds using dualBound_
4829  // for this we need clean basis so it is after factorize
4830  gutsOfSolution(NULL,NULL);
4831  numberFake_ =0; // Number of variables at fake bounds
4832  numberChanged_ =0; // Number of variables with changed costs
4833  changeBounds(true,NULL,objectiveChange);
4834
4835  problemStatus_ = -1;
4836  numberIterations_=0;
4837  factorization_->sparseThreshold(0);
4838  factorization_->goSparse();
4839
4840  int lastCleaned=0; // last time objective or bounds cleaned up
4841
4842  // number of times we have declared optimality
4843  numberTimesOptimal_=0;
4844
4845  // This says whether to restore things etc
4846  int factorType=0;
4847  /*
4848    Status of problem:
4849    0 - optimal
4850    1 - infeasible
4851    2 - unbounded
4852    -1 - iterating
4853    -2 - factorization wanted
4854    -3 - redo checking without factorization
4855    -4 - looks infeasible
4856
4857    BUT also from whileIterating return code is:
4858
4859   -1 iterations etc
4860   -2 inaccuracy
4861   -3 slight inaccuracy (and done iterations)
4862   +0 looks optimal (might be unbounded - but we will investigate)
4863   +1 looks infeasible
4864   +3 max iterations
4865
4866  */
4867
4868  int returnCode = 0;
4869
4870  int iRow,iColumn;
4871  while (problemStatus_<0) {
4872    // clear
4873    for (iRow=0;iRow<4;iRow++) {
4874      rowArray_[iRow]->clear();
4875    }   
4876   
4877    for (iColumn=0;iColumn<2;iColumn++) {
4878      columnArray_[iColumn]->clear();
4879    }   
4880
4881    // give matrix (and model costs and bounds a chance to be
4882    // refreshed (normally null)
4883    matrix_->refresh(this);
4884    // may factorize, checks if problem finished
4885    // should be able to speed this up on first time
4886    statusOfProblemInDual(lastCleaned,factorType,NULL,data,0);
4887
4888    // Say good factorization
4889    factorType=1;
4890
4891    // Do iterations
4892    if (problemStatus_<0) {
4893      double * givenPi=NULL;
4894      returnCode = whileIterating(givenPi,0);
4895      if ((!alwaysFinish&&returnCode<0)||returnCode==3) {
4896        if (returnCode!=3)
4897          assert (problemStatus_<0);
4898        returnCode=1;
4899        problemStatus_ = 3; 
4900        // can't say anything interesting - might as well return
4901#ifdef CLP_DEBUG
4902        printf("returning from fastDual after %d iterations with code %d\n",
4903               numberIterations_,returnCode);
4904#endif
4905        break;
4906      }
4907      returnCode=0;
4908    }
4909  }
4910
4911  // clear
4912  for (iRow=0;iRow<4;iRow++) {
4913    rowArray_[iRow]->clear();
4914  }   
4915 
4916  for (iColumn=0;iColumn<2;iColumn++) {
4917    columnArray_[iColumn]->clear();
4918  }   
4919  // Say not in fast dual
4920  specialOptions_ &= ~16384;
4921  assert(!numberFake_||((specialOptions_&(2048|4096))!=0&&dualBound_>1.0e8)
4922         ||returnCode||problemStatus_); // all bounds should be okay
4923  // Restore any saved stuff
4924  restoreData(data);
4925  dualBound_ = saveDualBound;
4926  return returnCode;
4927}
4928// This does first part of StrongBranching
4929ClpFactorization * 
4930ClpSimplexDual::setupForStrongBranching(char * arrays, int numberRows, int numberColumns)
4931{
4932  algorithm_ = -1;
4933  // put in standard form (and make row copy)
4934  // create modifiable copies of model rim and do optional scaling
4935  int startFinishOptions;
4936  /*  COIN_CLP_VETTED
4937      Looks safe for Cbc
4938  */
4939  bool safeOption = ((specialOptions_&COIN_CBC_USING_CLP)!=0);
4940  if((specialOptions_&4096)==0||(auxiliaryModel_&&!safeOption)) {
4941    startFinishOptions=0;
4942  } else {
4943    startFinishOptions=1+2+4;
4944  }
4945  createRim(7+8+16+32,true,startFinishOptions);
4946  // Do initial factorization
4947  // and set certain stuff
4948  // We can either set increasing rows so ...IsBasic gives pivot row
4949  // or we can just increment iBasic one by one
4950  // for now let ...iBasic give pivot row
4951  bool useFactorization=false;
4952  if ((startFinishOptions&2)!=0&&(whatsChanged_&(2+512))==2+512) {
4953    useFactorization=true; // Keep factorization if possible
4954    // switch off factorization if bad
4955    if (pivotVariable_[0]<0||factorization_->numberRows()!=numberRows_) 
4956      useFactorization=false;
4957  }
4958  if (!useFactorization) {
4959    factorization_->increasingRows(2);
4960    // row activities have negative sign
4961    factorization_->slackValue(-1.0);
4962    factorization_->zeroTolerance(1.0e-13);
4963
4964    int factorizationStatus = internalFactorize(0);
4965    if (factorizationStatus<0) {
4966      // some error
4967      // we should either debug or ignore
4968#ifndef NDEBUG
4969      printf("***** ClpDual strong branching factorization error - debug\n");
4970#endif
4971    } else if (factorizationStatus&&factorizationStatus<=numberRows_) {
4972      handler_->message(CLP_SINGULARITIES,messages_)
4973        <<factorizationStatus
4974        <<CoinMessageEol;
4975    }
4976  }
4977  double * arrayD = (double *) arrays;
4978  arrayD[0]=objectiveValue()*optimizationDirection_;
4979  double * saveSolution = arrayD+1;
4980  double * saveLower = saveSolution + (numberRows+numberColumns);
4981  double * saveUpper = saveLower + (numberRows+numberColumns);
4982  double * saveObjective = saveUpper + (numberRows+numberColumns);
4983  double * saveLowerOriginal = saveObjective + (numberRows+numberColumns);
4984  double * saveUpperOriginal = saveLowerOriginal + numberColumns;
4985  arrayD = saveUpperOriginal + numberColumns;
4986  int * savePivot = (int *) arrayD;
4987  int * whichRow = savePivot+numberRows;
4988  int * whichColumn = whichRow + 3*numberRows;
4989  int * arrayI = whichColumn + 2*numberColumns;
4990  unsigned char * saveStatus = (unsigned char *) (arrayI+1);
4991  // save stuff
4992  // save basis and solution
4993  CoinMemcpyN(solution_,
4994          numberRows_+numberColumns_,saveSolution);
4995  CoinMemcpyN(status_,numberColumns_+numberRows_,saveStatus);
4996  CoinMemcpyN(lower_,
4997          numberRows_+numberColumns_,saveLower);
4998  CoinMemcpyN(upper_,
4999          numberRows_+numberColumns_,saveUpper);
5000  CoinMemcpyN(cost_,
5001          numberRows_+numberColumns_,saveObjective);
5002  CoinMemcpyN(pivotVariable_, numberRows_,savePivot);
5003  return new ClpFactorization(*factorization_);
5004}
5005// This cleans up after strong branching
5006void 
5007ClpSimplexDual::cleanupAfterStrongBranching()
5008{
5009  int startFinishOptions;
5010  /*  COIN_CLP_VETTED
5011      Looks safe for Cbc
5012  */
5013  bool safeOption = ((specialOptions_&COIN_CBC_USING_CLP)!=0);
5014  if((specialOptions_&4096)==0||(auxiliaryModel_&&!safeOption)) {
5015    startFinishOptions=0;
5016  } else {
5017    startFinishOptions=1+2+4;
5018  }
5019  if ((startFinishOptions&1)==0) {
5020    deleteRim(1);
5021    whatsChanged_=0;
5022  } else {
5023    // Original factorization will have been put back by last loop
5024    //delete factorization_;
5025    //factorization_ = new ClpFactorization(saveFactorization);
5026    deleteRim(0);
5027    // mark all as current
5028    whatsChanged_ = 0xffff;
5029  }
5030}
5031/* Checks number of variables at fake bounds.  This is used by fastDual
5032   so can exit gracefully before end */
5033int 
5034ClpSimplexDual::numberAtFakeBound()
5035{
5036  int iSequence;
5037  int numberFake=0;
5038 
5039  for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
5040    FakeBound bound = getFakeBound(iSequence);
5041    switch(getStatus(iSequence)) {
5042
5043    case basic:
5044      break;
5045    case isFree:
5046    case superBasic:
5047    case ClpSimplex::isFixed:
5048      assert (bound==noFake);
5049      break;
5050    case atUpperBound:
5051      if (bound==upperFake||bound==bothFake)
5052        numberFake++;
5053      break;
5054    case atLowerBound:
5055      if (bound==lowerFake||bound==bothFake)
5056        numberFake++;
5057      break;
5058    }
5059  }
5060  numberFake_ = numberFake;
5061  return numberFake;
5062}
5063/* Pivot out a variable and choose an incoing one.  Assumes dual
5064   feasible - will not go through a reduced cost. 
5065   Returns step length in theta
5066   Returns ray in ray_ (or NULL if no pivot)
5067   Return codes as before but -1 means no acceptable pivot
5068*/
5069int 
5070ClpSimplexDual::pivotResult()
5071{
5072  abort();
5073  return 0;
5074}
5075/*
5076   Row array has row part of pivot row
5077   Column array has column part.
5078   This is used in dual values pass
5079*/
5080void
5081ClpSimplexDual::checkPossibleValuesMove(CoinIndexedVector * rowArray,
5082                                        CoinIndexedVector * columnArray,
5083                                        double acceptablePivot)
5084{
5085  double * work;
5086  int number;
5087  int * which;
5088  int iSection;
5089
5090  double tolerance = dualTolerance_*1.001;
5091
5092  double thetaDown = 1.0e31;
5093  double changeDown ;
5094  double thetaUp = 1.0e31;
5095  double bestAlphaDown = acceptablePivot*0.99999;
5096  double bestAlphaUp = acceptablePivot*0.99999;
5097  int sequenceDown =-1;
5098  int sequenceUp = sequenceOut_;
5099
5100  double djBasic = dj_[sequenceOut_];
5101  if (djBasic>0.0) {
5102    // basic at lower bound so directionOut_ 1 and -1 in pivot row
5103    // dj will go to zero on other way
5104    thetaUp = djBasic;
5105    changeDown = -lower_[sequenceOut_];
5106  } else {
5107    // basic at upper bound so directionOut_ -1 and 1 in pivot row
5108    // dj will go to zero on other way
5109    thetaUp = -djBasic;
5110    changeDown = upper_[sequenceOut_];
5111  }
5112  bestAlphaUp = 1.0;
5113  int addSequence;
5114
5115  double alphaUp=0.0;
5116  double alphaDown=0.0;
5117
5118  for (iSection=0;iSection<2;iSection++) {
5119
5120    int i;
5121    if (!iSection) {
5122      work = rowArray->denseVector();
5123      number = rowArray->getNumElements();
5124      which = rowArray->getIndices();
5125      addSequence = numberColumns_;
5126    } else {
5127      work = columnArray->denseVector();
5128      number = columnArray->getNumElements();
5129      which = columnArray->getIndices();
5130      addSequence = 0;
5131    }
5132   
5133    for (i=0;i<number;i++) {
5134      int iSequence = which[i];
5135      int iSequence2 = iSequence + addSequence;
5136      double alpha;
5137      double oldValue;
5138      double value;
5139
5140      switch(getStatus(iSequence2)) {
5141         
5142      case basic: 
5143        break;
5144      case ClpSimplex::isFixed:
5145        alpha = work[i];
5146        changeDown += alpha*upper_[iSequence2];
5147        break;
5148      case isFree:
5149      case superBasic:
5150        alpha = work[i];
5151        // dj must be effectively zero as dual feasible
5152        if (fabs(alpha)>bestAlphaUp) {
5153          thetaDown = 0.0;
5154          thetaUp = 0.0;
5155          bestAlphaDown = fabs(alpha);
5156          bestAlphaUp = bestAlphaUp;
5157          sequenceDown =iSequence2;
5158          sequenceUp = sequenceDown;
5159          alphaUp = alpha;
5160          alphaDown = alpha;
5161        }
5162        break;
5163      case atUpperBound:
5164        alpha = work[i];
5165        oldValue = dj_[iSequence2];
5166        changeDown += alpha*upper_[iSequence2];
5167        if (alpha>=acceptablePivot) {
5168          // might do other way
5169          value = oldValue+thetaUp*alpha;
5170          if (value>-tolerance) {
5171            if (value>tolerance||fabs(alpha)>bestAlphaUp) {
5172              thetaUp = -oldValue/alpha;
5173              bestAlphaUp = fabs(alpha);
5174              sequenceUp = iSequence2;
5175              alphaUp = alpha;
5176            }
5177          }
5178        } else if (alpha<=-acceptablePivot) {
5179          // might do this way
5180          value = oldValue-thetaDown*alpha;
5181          if (value>-tolerance) {
5182            if (value>tolerance||fabs(alpha)>bestAlphaDown) {
5183              thetaDown = oldValue/alpha;
5184              bestAlphaDown = fabs(alpha);
5185              sequenceDown = iSequence2;
5186              alphaDown = alpha;
5187            }
5188          }
5189        }
5190        break;
5191      case atLowerBound:
5192        alpha = work[i];
5193        oldValue = dj_[iSequence2];
5194        changeDown += alpha*lower_[iSequence2];
5195        if (alpha<=-acceptablePivot) {
5196          // might do other way
5197          value = oldValue+thetaUp*alpha;
5198          if (value<tolerance) {
5199            if (value<-tolerance||fabs(alpha)>bestAlphaUp) {
5200              thetaUp = -oldValue/alpha;
5201              bestAlphaUp = fabs(alpha);
5202              sequenceUp = iSequence2;
5203              alphaUp = alpha;
5204            }
5205          }
5206        } else if (alpha>=acceptablePivot) {
5207          // might do this way
5208          value = oldValue-thetaDown*alpha;
5209          if (value<tolerance) {
5210            if (value<-tolerance||fabs(alpha)>bestAlphaDown) {
5211              thetaDown = oldValue/alpha;
5212              bestAlphaDown = fabs(alpha);
5213              sequenceDown = iSequence2;
5214              alphaDown = alpha;
5215            }
5216          }
5217        }
5218        break;
5219      }
5220    }
5221  }
5222  thetaUp *= -1.0;
5223  double changeUp = -thetaUp * changeDown;
5224  changeDown = -thetaDown*changeDown;
5225  if (CoinMax(fabs(thetaDown),fabs(thetaUp))<1.0e-8) {
5226    // largest
5227    if (fabs(alphaDown)<fabs(alphaUp)) {
5228      sequenceDown =-1;
5229    }
5230  }
5231  // choose
5232  sequenceIn_=-1;
5233  if (changeDown>changeUp&&sequenceDown>=0) {
5234    theta_ = thetaDown;
5235    if (fabs(changeDown)<1.0e30)
5236      sequenceIn_ = sequenceDown;
5237    alpha_ = alphaDown;
5238#ifdef CLP_DEBUG
5239    if ((handler_->logLevel()&32))
5240      printf("predicted way - dirout %d, change %g,%g theta %g\n",
5241             directionOut_,changeDown,changeUp,theta_);
5242#endif
5243  } else {
5244    theta_ = thetaUp;
5245    if (fabs(changeUp)<1.0e30)
5246      sequenceIn_ = sequenceUp;
5247    alpha_ = alphaUp;
5248    if (sequenceIn_!=sequenceOut_) {
5249#ifdef CLP_DEBUG
5250      if ((handler_->logLevel()&32))
5251        printf("opposite way - dirout %d, change %g,%g theta %g\n",
5252               directionOut_,changeDown,changeUp,theta_);
5253#endif
5254    } else {
5255#ifdef CLP_DEBUG
5256      if ((handler_->logLevel()&32))
5257        printf("opposite way to zero dj - dirout %d, change %g,%g theta %g\n",
5258               directionOut_,changeDown,changeUp,theta_);
5259#endif
5260    }
5261  }
5262  if (sequenceIn_>=0) {
5263    lowerIn_ = lower_[sequenceIn_];
5264    upperIn_ = upper_[sequenceIn_];
5265    valueIn_ = solution_[sequenceIn_];
5266    dualIn_ = dj_[sequenceIn_];
5267
5268    if (alpha_<0.0) {
5269      // as if from upper bound
5270      directionIn_=-1;
5271      upperIn_=valueIn_;
5272    } else {
5273      // as if from lower bound
5274      directionIn_=1;
5275      lowerIn_=valueIn_;
5276    }
5277  }
5278}
5279/*
5280   Row array has row part of pivot row
5281   Column array has column part.
5282   This is used in cleanup
5283*/
5284void
5285ClpSimplexDual::checkPossibleCleanup(CoinIndexedVector * rowArray,
5286                                        CoinIndexedVector * columnArray,
5287                                        double acceptablePivot)
5288{
5289  double * work;
5290  int number;
5291  int * which;
5292  int iSection;
5293
5294  double tolerance = dualTolerance_*1.001;
5295
5296  double thetaDown = 1.0e31;
5297  double thetaUp = 1.0e31;
5298  double bestAlphaDown = acceptablePivot*10.0;
5299  double bestAlphaUp = acceptablePivot*10.0;
5300  int sequenceDown =-1;
5301  int sequenceUp = -1;
5302
5303  double djSlack = dj_[pivotRow_];
5304  if (getRowStatus(pivotRow_)==basic)
5305    djSlack=COIN_DBL_MAX;
5306  if (fabs(djSlack)<tolerance)
5307    djSlack=0.0;
5308  int addSequence;
5309
5310  double alphaUp=0.0;
5311  double alphaDown=0.0;
5312  for (iSection=0;iSection<2;iSection++) {
5313
5314    int i;
5315    if (!iSection) {
5316      work = rowArray->denseVector();
5317      number = rowArray->getNumElements();
5318      which = rowArray->getIndices();
5319      addSequence = numberColumns_;
5320    } else {
5321      work = columnArray->denseVector();
5322      number = columnArray->getNumElements();
5323      which = columnArray->getIndices();
5324      addSequence = 0;
5325    }
5326   
5327    for (i=0;i<number;i++) {
5328      int iSequence = which[i];
5329      int iSequence2 = iSequence + addSequence;
5330      double alpha;
5331      double oldValue;
5332      double value;
5333
5334      switch(getStatus(iSequence2)) {
5335         
5336      case basic: 
5337        break;
5338      case ClpSimplex::isFixed:
5339        alpha = work[i];
5340        if (addSequence) {
5341          printf("possible - pivot row %d this %d\n",pivotRow_,iSequence);
5342          oldValue = dj_[iSequence2];
5343          if (alpha<=-acceptablePivot) {
5344            // might do other way
5345            value = oldValue+thetaUp*alpha;
5346            if (value<tolerance) {
5347              if (value<-tolerance||fabs(alpha)>bestAlphaUp) {
5348                thetaUp = -oldValue/alpha;
5349                bestAlphaUp = fabs(alpha);
5350                sequenceUp = iSequence2;
5351                alphaUp = alpha;
5352              }
5353            }
5354          } else if (alpha>=acceptablePivot) {
5355            // might do this way
5356            value = oldValue-thetaDown*alpha;
5357            if (value<tolerance) {
5358              if (value<-tolerance||fabs(alpha)>bestAlphaDown) {
5359                thetaDown = oldValue/alpha;
5360                bestAlphaDown = fabs(alpha);
5361                sequenceDown = iSequence2;
5362                alphaDown = alpha;
5363              }
5364            }
5365          }
5366        }
5367        break;
5368      case isFree:
5369      case superBasic:
5370        alpha = work[i];
5371        // dj must be effectively zero as dual feasible
5372        if (fabs(alpha)>bestAlphaUp) {
5373          thetaDown = 0.0;
5374          thetaUp = 0.0;
5375          bestAlphaDown = fabs(alpha);
5376          bestAlphaUp = bestAlphaUp;
5377          sequenceDown =iSequence2;
5378          sequenceUp = sequenceDown;
5379          alphaUp = alpha;
5380          alphaDown = alpha;
5381        }
5382        break;
5383      case atUpperBound:
5384        alpha = work[i];
5385        oldValue = dj_[iSequence2];
5386        if (alpha>=acceptablePivot) {
5387          // might do other way
5388          value = oldValue+thetaUp*alpha;
5389          if (value>-tolerance) {
5390            if (value>tolerance||fabs(alpha)>bestAlphaUp) {
5391              thetaUp = -oldValue/alpha;
5392              bestAlphaUp = fabs(alpha);
5393              sequenceUp = iSequence2;
5394              alphaUp = alpha;
5395            }
5396          }
5397        } else if (alpha<=-acceptablePivot) {
5398          // might do this way
5399          value = oldValue-thetaDown*alpha;
5400          if (value>-tolerance) {
5401            if (value>tolerance||fabs(alpha)>bestAlphaDown) {
5402              thetaDown = oldValue/alpha;
5403              bestAlphaDown = fabs(alpha);
5404              sequenceDown = iSequence2;
5405              alphaDown = alpha;
5406            }
5407          }
5408        }
5409        break;
5410      case atLowerBound:
5411        alpha = work[i];
5412        oldValue = dj_[iSequence2];
5413        if (alpha<=-acceptablePivot) {
5414          // might do other way
5415          value = oldValue+thetaUp*alpha;
5416          if (value<tolerance) {
5417            if (value<-tolerance||fabs(alpha)>bestAlphaUp) {
5418              thetaUp = -oldValue/alpha;
5419              bestAlphaUp = fabs(alpha);
5420              sequenceUp = iSequence2;
5421              alphaUp = alpha;
5422            }
5423          }
5424        } else if (alpha>=acceptablePivot) {
5425          // might do this way
5426          value = oldValue-thetaDown*alpha;
5427          if (value<tolerance) {
5428            if (value<-tolerance||fabs(alpha)>bestAlphaDown) {
5429              thetaDown = oldValue/alpha;
5430              bestAlphaDown = fabs(alpha);
5431              sequenceDown = iSequence2;
5432              alphaDown = alpha;
5433            }
5434          }
5435        }
5436        break;
5437      }
5438    }
5439  }
5440  thetaUp *= -1.0;
5441  // largest
5442  if (bestAlphaDown<bestAlphaUp) 
5443    sequenceDown =-1;
5444  else
5445    sequenceUp=-1;
5446
5447  sequenceIn_=-1;
5448 
5449  if (sequenceDown>=0) {
5450    theta_ = thetaDown;
5451    sequenceIn_ = sequenceDown;
5452    alpha_ = alphaDown;
5453#ifdef CLP_DEBUG
5454    if ((handler_->logLevel()&32))
5455      printf("predicted way - dirout %d, theta %g\n",
5456             directionOut_,theta_);
5457#endif
5458  } else if (sequenceUp>=0) {
5459    theta_ = thetaUp;
5460    sequenceIn_ = sequenceUp;
5461    alpha_ = alphaUp;
5462#ifdef CLP_DEBUG
5463    if ((handler_->logLevel()&32))
5464      printf("opposite way - dirout %d,theta %g\n",
5465             directionOut_,theta_);
5466#endif
5467  }
5468  if (sequenceIn_>=0) {
5469    lowerIn_ = lower_[sequenceIn_];
5470    upperIn_ = upper_[sequenceIn_];
5471    valueIn_ = solution_[sequenceIn_];
5472    dualIn_ = dj_[sequenceIn_];
5473
5474    if (alpha_<0.0) {
5475      // as if from upper bound
5476      directionIn_=-1;
5477      upperIn_=valueIn_;
5478    } else {
5479      // as if from lower bound
5480      directionIn_=1;
5481      lowerIn_=valueIn_;
5482    }
5483  }
5484}
5485/*
5486   This sees if we can move duals in dual values pass.
5487   This is done before any pivoting
5488*/
5489void ClpSimplexDual::doEasyOnesInValuesPass(double * dj)
5490{
5491  // Get column copy
5492  CoinPackedMatrix * columnCopy = matrix();
5493  // Get a row copy in standard format
5494  CoinPackedMatrix copy;
5495  copy.reverseOrderedCopyOf(*columnCopy);
5496  // get matrix data pointers
5497  const int * column = copy.getIndices();
5498  const CoinBigIndex * rowStart = copy.getVectorStarts();
5499  const int * rowLength = copy.getVectorLengths(); 
5500  const double * elementByRow = copy.getElements();
5501  double tolerance = dualTolerance_*1.001;
5502
5503  int iRow;
5504#ifdef CLP_DEBUG
5505  {
5506    double value5=0.0;
5507    int i;
5508    for (i=0;i<numberRows_+numberColumns_;i++) {
5509      if (dj[i]<-1.0e-6)
5510        value5 += dj[i]*upper_[i];
5511      else if (dj[i] >1.0e-6)
5512        value5 += dj[i]*lower_[i];
5513    }
5514    printf("Values objective Value before %g\n",value5);
5515  }
5516#endif
5517  // for scaled row
5518  double * scaled=NULL;
5519  if (rowScale_)
5520    scaled = new double[numberColumns_];
5521  for (iRow=0;iRow<numberRows_;iRow++) {
5522
5523    int iSequence = iRow + numberColumns_;
5524    double djBasic = dj[iSequence];
5525    if (getRowStatus(iRow)==basic&&fabs(djBasic)>tolerance) {
5526
5527      double changeUp ;
5528      // always -1 in pivot row
5529      if (djBasic>0.0) {
5530        // basic at lower bound
5531        changeUp = -lower_[iSequence];
5532      } else {
5533        // basic at upper bound
5534        changeUp = upper_[iSequence];
5535      }
5536      bool canMove=true;
5537      int i;
5538      const double * thisElements = elementByRow + rowStart[iRow]; 
5539      const int * thisIndices = column+rowStart[iRow];
5540      if (rowScale_) {
5541        assert (!auxiliaryModel_);
5542        // scale row
5543        double scale = rowScale_[iRow];
5544        for (i = 0;i<rowLength[iRow];i++) {
5545          int iColumn = thisIndices[i];
5546          double alpha = thisElements[i];
5547          scaled[i] = scale*alpha*columnScale_[iColumn];
5548        }
5549        thisElements = scaled;
5550      }
5551      for (i = 0;i<rowLength[iRow];i++) {
5552        int iColumn = thisIndices[i];
5553        double alpha = thisElements[i];
5554        double oldValue = dj[iColumn];;
5555        double value;
5556
5557        switch(getStatus(iColumn)) {
5558         
5559        case basic:
5560          if (dj[iColumn]<-tolerance&&
5561              fabs(solution_[iColumn]-upper_[iColumn])<1.0e-8) {
5562            // at ub
5563            changeUp += alpha*upper_[iColumn];
5564            // might do other way
5565            value = oldValue+djBasic*alpha;
5566            if (value>tolerance) 
5567              canMove=false;
5568          } else if (dj[iColumn]>tolerance&&
5569              fabs(solution_[iColumn]-lower_[iColumn])<1.0e-8) {
5570            changeUp += alpha*lower_[iColumn];
5571            // might do other way
5572            value = oldValue+djBasic*alpha;
5573            if (value<-tolerance)
5574              canMove=false;
5575          } else {
5576            canMove=false;
5577          }
5578          break;
5579        case ClpSimplex::isFixed:
5580          changeUp += alpha*upper_[iColumn];
5581          break;
5582        case isFree:
5583        case superBasic:
5584          canMove=false;
5585        break;
5586      case atUpperBound:
5587        changeUp += alpha*upper_[iColumn];
5588        // might do other way
5589        value = oldValue+djBasic*alpha;
5590        if (value>tolerance) 
5591          canMove=false;
5592        break;
5593      case atLowerBound:
5594        changeUp += alpha*lower_[iColumn];
5595        // might do other way
5596        value = oldValue+djBasic*alpha;
5597        if (value<-tolerance)
5598          canMove=false;
5599        break;
5600        }
5601      }
5602      if (canMove) {
5603        if (changeUp*djBasic>1.0e-12||fabs(changeUp)<1.0e-8) {
5604          // move
5605          for (i = 0;i<rowLength[iRow];i++) {
5606            int iColumn = thisIndices[i];
5607            double alpha = thisElements[i];
5608            dj[iColumn] += djBasic * alpha;
5609          }
5610          dj[iSequence]=0.0;
5611#ifdef CLP_DEBUG
5612          {
5613            double value5=0.0;
5614            int i;
5615            for (i=0;i<numberRows_+numberColumns_;i++) {
5616              if (dj[i]<-1.0e-6)
5617                value5 += dj[i]*upper_[i];
5618              else if (dj[i] >1.0e-6)
5619                value5 += dj[i]*lower_[i];
5620            }
5621            printf("Values objective Value after row %d old dj %g %g\n",
5622                   iRow,djBasic,value5);
5623          }
5624#endif
5625        }
5626      }
5627    }     
5628  }
5629  delete [] scaled;
5630}
5631int
5632ClpSimplexDual::nextSuperBasic()
5633{
5634  if (firstFree_>=0) {
5635    int returnValue=firstFree_;
5636    int iColumn=firstFree_+1;
5637    for (;iColumn<numberRows_+numberColumns_;iColumn++) {
5638      if (getStatus(iColumn)==isFree) 
5639        if (fabs(dj_[iColumn])>1.0e2*dualTolerance_) 
5640          break;
5641    }
5642    firstFree_ = iColumn;
5643    if (firstFree_==numberRows_+numberColumns_)
5644      firstFree_=-1;
5645    return returnValue;
5646  } else {
5647    return -1;
5648  }
5649}
Note: See TracBrowser for help on using the repository browser.