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

Last change on this file since 1287 was 1287, checked in by forrest, 11 years ago

changes for cbc event handler and multiple factorizations

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