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

Last change on this file since 1429 was 1429, checked in by forrest, 10 years ago

fixes etc

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