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

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

left a printf in

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