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

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

changes to try and make faster

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