source: stable/1.11/Clp/src/ClpSimplexDual.cpp @ 1507

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

fix for free variables

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