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

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

add ids

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