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

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

fix for ray and warning messages

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