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

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

fix assert on scaling

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