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

Last change on this file since 1301 was 1301, checked in by forrest, 13 years ago

try and fix fake feasible for Scip

  • 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                printf("XXX Infeas ? %d inf summing to %g\n",numberPrimalInfeasibilities_,
1826                       sumPrimalInfeasibilities_);
1827                problemStatus_=-1;
1828                returnCode=-2;
1829              }
1830#ifndef NDEBUG
1831              memcpy(solution_,comp,nTotal*sizeof(double));
1832              delete [] comp;
1833#endif
1834            }
1835#endif
1836            if (!problemStatus_) {
1837              // make it look OK
1838              numberPrimalInfeasibilities_=0;
1839              sumPrimalInfeasibilities_=0.0;
1840              numberDualInfeasibilities_=0;
1841              sumDualInfeasibilities_=0.0;
1842              // May be perturbed
1843              if (perturbation_==101||numberChanged_) {
1844                perturbation_=102; // stop any perturbations
1845                //double changeCost;
1846                //changeBounds(1,NULL,changeCost);
1847                createRim4(false);
1848                // make sure duals are current
1849                computeDuals(givenDuals);
1850                checkDualSolution();
1851                if (numberDualInfeasibilities_) {
1852                  problemStatus_=10; // was -3;
1853                } else {
1854                  computeObjectiveValue(true);
1855                }
1856              } else if (numberPivots) {
1857                computeObjectiveValue(true);
1858              } 
1859              if (numberPivots<-1000) {
1860                // objective may be wrong
1861                objectiveValue_ = innerProduct(cost_,numberColumns_+numberRows_,solution_);
1862                objectiveValue_ += objective_->nonlinearOffset();
1863                objectiveValue_ /= (objectiveScale_*rhsScale_);
1864                if ((specialOptions_&16384)==0) {
1865                  // and dual_ may be wrong (i.e. for fixed or basic)
1866                  CoinIndexedVector * arrayVector = rowArray_[1];
1867                  arrayVector->clear();
1868                  int iRow;
1869                  double * array = arrayVector->denseVector();
1870                  /* Use dual_ instead of array
1871                     Even though dual_ is only numberRows_ long this is
1872                     okay as gets permuted to longer rowArray_[2]
1873                  */
1874                  arrayVector->setDenseVector(dual_);
1875                  int * index = arrayVector->getIndices();
1876                  int number=0;
1877                  for (iRow=0;iRow<numberRows_;iRow++) {
1878                    int iPivot=pivotVariable_[iRow];
1879                    double value = cost_[iPivot];
1880                    dual_[iRow]=value;
1881                    if (value) {
1882                      index[number++]=iRow;
1883                    }
1884                  }
1885                  arrayVector->setNumElements(number);
1886                  // Extended duals before "updateTranspose"
1887                  matrix_->dualExpanded(this,arrayVector,NULL,0);
1888                  // Btran basic costs
1889                  rowArray_[2]->clear();
1890                  factorization_->updateColumnTranspose(rowArray_[2],arrayVector);
1891                  // and return vector
1892                  arrayVector->setDenseVector(array);
1893                }
1894              }
1895              sumPrimalInfeasibilities_=0.0;
1896            }
1897            if ((specialOptions_&(1024+16384))!=0&&!problemStatus_) {
1898              CoinIndexedVector * arrayVector = rowArray_[1];
1899              arrayVector->clear();
1900              double * rhs = arrayVector->denseVector();
1901              times(1.0,solution_,rhs);
1902#ifdef CHECK_ACCURACY
1903              bool bad=false;
1904#endif
1905              bool bad2=false;
1906              int i;
1907              for ( i=0;i<numberRows_;i++) {
1908                if (rhs[i]<rowLowerWork_[i]-primalTolerance_||
1909                    rhs[i]>rowUpperWork_[i]+primalTolerance_) {
1910                  bad2=true;
1911#ifdef CHECK_ACCURACY
1912                  printf("row %d out of bounds %g, %g correct %g bad %g\n",i,
1913                         rowLowerWork_[i],rowUpperWork_[i],
1914                         rhs[i],rowActivityWork_[i]);
1915#endif
1916                } else if (fabs(rhs[i]-rowActivityWork_[i])>1.0e-3) {
1917#ifdef CHECK_ACCURACY
1918                  bad=true;
1919                  printf("row %d correct %g bad %g\n",i,rhs[i],rowActivityWork_[i]);
1920#endif
1921                }
1922                rhs[i]=0.0;
1923              }
1924              for ( i=0;i<numberColumns_;i++) {
1925                if (solution_[i]<columnLowerWork_[i]-primalTolerance_||
1926                    solution_[i]>columnUpperWork_[i]+primalTolerance_) {
1927                  bad2=true;
1928#ifdef CHECK_ACCURACY
1929                  printf("column %d out of bounds %g, %g correct %g bad %g\n",i,
1930                         columnLowerWork_[i],columnUpperWork_[i],
1931                         solution_[i],columnActivityWork_[i]);
1932#endif
1933                }
1934              }
1935              if (bad2) {
1936                problemStatus_=-3;
1937                returnCode=-2;
1938                // Force to re-factorize early next time
1939                int numberPivots = factorization_->pivots();
1940                forceFactorization_=CoinMin(forceFactorization_,(numberPivots+1)>>1);
1941              }
1942            }
1943          }
1944        }
1945      } else {
1946        problemStatus_=-3;
1947        returnCode=-2;
1948        // Force to re-factorize early next time
1949        int numberPivots = factorization_->pivots();
1950        forceFactorization_=CoinMin(forceFactorization_,(numberPivots+1)>>1);
1951      }
1952      break;
1953    }
1954  }
1955  if (givenDuals) {
1956    CoinMemcpyN(dj_,numberRows_+numberColumns_,givenDuals);
1957    // get rid of any values pass array
1958    delete [] candidateList;
1959  }
1960  delete [] dubiousWeights;
1961  return returnCode;
1962}
1963/* The duals are updated by the given arrays.
1964   Returns number of infeasibilities.
1965   rowArray and columnarray will have flipped
1966   The output vector has movement (row length array) */
1967int
1968ClpSimplexDual::updateDualsInDual(CoinIndexedVector * rowArray,
1969                                  CoinIndexedVector * columnArray,
1970                                  CoinIndexedVector * outputArray,
1971                                  double theta,
1972                                  double & objectiveChange,
1973                                  bool fullRecompute)
1974{
1975 
1976  outputArray->clear();
1977 
1978 
1979  int numberInfeasibilities=0;
1980  int numberRowInfeasibilities=0;
1981 
1982  // get a tolerance
1983  double tolerance = dualTolerance_;
1984  // we can't really trust infeasibilities if there is dual error
1985  double error = CoinMin(1.0e-2,largestDualError_);
1986  // allow tolerance at least slightly bigger than standard
1987  tolerance = tolerance +  error;
1988 
1989  double changeObj=0.0;
1990
1991  // Coding is very similar but we can save a bit by splitting
1992  // Do rows
1993  if (!fullRecompute) {
1994    int i;
1995    double * reducedCost = djRegion(0);
1996    const double * lower = lowerRegion(0);
1997    const double * upper = upperRegion(0);
1998    const double * cost = costRegion(0);
1999    double * work;
2000    int number;
2001    int * which;
2002    assert(rowArray->packedMode());
2003    work = rowArray->denseVector();
2004    number = rowArray->getNumElements();
2005    which = rowArray->getIndices();
2006    for (i=0;i<number;i++) {
2007      int iSequence = which[i];
2008      double alphaI = work[i];
2009      work[i]=0.0;
2010     
2011      Status status = getStatus(iSequence+numberColumns_);
2012      // more likely to be at upper bound ?
2013      if (status==atUpperBound) {
2014
2015        double value = reducedCost[iSequence]-theta*alphaI;
2016        reducedCost[iSequence]=value;
2017
2018        if (value>tolerance) {
2019          // to lower bound (if swap)
2020          which[numberInfeasibilities++]=iSequence;
2021          double movement = lower[iSequence]-upper[iSequence];
2022          assert (fabs(movement)<1.0e30);
2023#ifdef CLP_DEBUG
2024          if ((handler_->logLevel()&32))
2025            printf("%d %d, new dj %g, alpha %g, movement %g\n",
2026                   0,iSequence,value,alphaI,movement);
2027#endif
2028          changeObj += movement*cost[iSequence];
2029          outputArray->quickAdd(iSequence,-movement);
2030        }
2031      } else if (status==atLowerBound) {
2032
2033        double value = reducedCost[iSequence]-theta*alphaI;
2034        reducedCost[iSequence]=value;
2035
2036        if (value<-tolerance) {
2037          // to upper bound
2038          which[numberInfeasibilities++]=iSequence;
2039          double movement = upper[iSequence] - lower[iSequence];
2040          assert (fabs(movement)<1.0e30);
2041#ifdef CLP_DEBUG
2042          if ((handler_->logLevel()&32))
2043            printf("%d %d, new dj %g, alpha %g, movement %g\n",
2044                   0,iSequence,value,alphaI,movement);
2045#endif
2046          changeObj += movement*cost[iSequence];
2047          outputArray->quickAdd(iSequence,-movement);
2048        }
2049      }
2050    }
2051  } else  {
2052    double * solution = solutionRegion(0);
2053    double * reducedCost = djRegion(0);
2054    const double * lower = lowerRegion(0);
2055    const double * upper = upperRegion(0);
2056    const double * cost = costRegion(0);
2057    int * which;
2058    which = rowArray->getIndices();
2059    int iSequence;
2060    for (iSequence=0;iSequence<numberRows_;iSequence++) {
2061      double value = reducedCost[iSequence];
2062     
2063      Status status = getStatus(iSequence+numberColumns_);
2064      // more likely to be at upper bound ?
2065      if (status==atUpperBound) {
2066        double movement=0.0;
2067
2068        if (value>tolerance) {
2069          // to lower bound (if swap)
2070          // put back alpha
2071          which[numberInfeasibilities++]=iSequence;
2072          movement = lower[iSequence]-upper[iSequence];
2073          changeObj += movement*cost[iSequence];
2074          outputArray->quickAdd(iSequence,-movement);
2075          assert (fabs(movement)<1.0e30);
2076        } else if (value>-tolerance) {
2077          // at correct bound but may swap
2078          FakeBound bound = getFakeBound(iSequence+numberColumns_);
2079          if (bound==ClpSimplexDual::upperFake) {
2080            movement = lower[iSequence]-upper[iSequence];
2081            assert (fabs(movement)<1.0e30);
2082            setStatus(iSequence+numberColumns_,atLowerBound);
2083            solution[iSequence] = lower[iSequence];
2084            changeObj += movement*cost[iSequence];
2085            numberFake_--;
2086            setFakeBound(iSequence+numberColumns_,noFake);
2087          }
2088        }
2089      } else if (status==atLowerBound) {
2090        double movement=0.0;
2091
2092        if (value<-tolerance) {
2093          // to upper bound
2094          // put back alpha
2095          which[numberInfeasibilities++]=iSequence;
2096          movement = upper[iSequence] - lower[iSequence];
2097          assert (fabs(movement)<1.0e30);
2098          changeObj += movement*cost[iSequence];
2099          outputArray->quickAdd(iSequence,-movement);
2100        } else if (value<tolerance) {
2101          // at correct bound but may swap
2102          FakeBound bound = getFakeBound(iSequence+numberColumns_);
2103          if (bound==ClpSimplexDual::lowerFake) {
2104            movement = upper[iSequence]-lower[iSequence];
2105            assert (fabs(movement)<1.0e30);
2106            setStatus(iSequence+numberColumns_,atUpperBound);
2107            solution[iSequence] = upper[iSequence];
2108            changeObj += movement*cost[iSequence];
2109            numberFake_--;
2110            setFakeBound(iSequence+numberColumns_,noFake);
2111          }
2112        }
2113      }
2114    }
2115  }
2116
2117  // Do columns
2118  if (!fullRecompute) {
2119    int i;
2120    double * reducedCost = djRegion(1);
2121    const double * lower = lowerRegion(1);
2122    const double * upper = upperRegion(1);
2123    const double * cost = costRegion(1);
2124    double * work;
2125    int number;
2126    int * which;
2127    // set number of infeasibilities in row array
2128    numberRowInfeasibilities=numberInfeasibilities;
2129    rowArray->setNumElements(numberInfeasibilities);
2130    numberInfeasibilities=0;
2131    work = columnArray->denseVector();
2132    number = columnArray->getNumElements();
2133    which = columnArray->getIndices();
2134   
2135    for (i=0;i<number;i++) {
2136      int iSequence = which[i];
2137      double alphaI = work[i];
2138      work[i]=0.0;
2139     
2140      Status status = getStatus(iSequence);
2141      if (status==atLowerBound) {
2142        double value = reducedCost[iSequence]-theta*alphaI;
2143        reducedCost[iSequence]=value;
2144        double movement=0.0;
2145
2146        if (value<-tolerance) {
2147          // to upper bound
2148          which[numberInfeasibilities++]=iSequence;
2149          movement = upper[iSequence] - lower[iSequence];
2150          assert (fabs(movement)<1.0e30);
2151#ifdef CLP_DEBUG
2152          if ((handler_->logLevel()&32))
2153            printf("%d %d, new dj %g, alpha %g, movement %g\n",
2154                   1,iSequence,value,alphaI,movement);
2155#endif
2156          changeObj += movement*cost[iSequence];
2157          matrix_->add(this,outputArray,iSequence,movement);
2158        }
2159      } else if (status==atUpperBound) {
2160        double value = reducedCost[iSequence]-theta*alphaI;
2161        reducedCost[iSequence]=value;
2162        double movement=0.0;
2163
2164        if (value>tolerance) {
2165          // to lower bound (if swap)
2166          which[numberInfeasibilities++]=iSequence;
2167          movement = lower[iSequence]-upper[iSequence];
2168          assert (fabs(movement)<1.0e30);
2169#ifdef CLP_DEBUG
2170          if ((handler_->logLevel()&32))
2171            printf("%d %d, new dj %g, alpha %g, movement %g\n",
2172                   1,iSequence,value,alphaI,movement);
2173#endif
2174          changeObj += movement*cost[iSequence];
2175          matrix_->add(this,outputArray,iSequence,movement);
2176        }
2177      } else if (status==isFree) {
2178        double value = reducedCost[iSequence]-theta*alphaI;
2179        reducedCost[iSequence]=value;
2180      }
2181    }
2182  } else {
2183    double * solution = solutionRegion(1);
2184    double * reducedCost = djRegion(1);
2185    const double * lower = lowerRegion(1);
2186    const double * upper = upperRegion(1);
2187    const double * cost = costRegion(1);
2188    int * which;
2189    // set number of infeasibilities in row array
2190    numberRowInfeasibilities=numberInfeasibilities;
2191    rowArray->setNumElements(numberInfeasibilities);
2192    numberInfeasibilities=0;
2193    which = columnArray->getIndices();
2194    int iSequence;
2195    for (iSequence=0;iSequence<numberColumns_;iSequence++) {
2196      double value = reducedCost[iSequence];
2197     
2198      Status status = getStatus(iSequence);
2199      if (status==atLowerBound) {
2200        double movement=0.0;
2201
2202        if (value<-tolerance) {
2203          // to upper bound
2204          // put back alpha
2205          which[numberInfeasibilities++]=iSequence;
2206          movement = upper[iSequence] - lower[iSequence];
2207          assert (fabs(movement)<1.0e30);
2208          changeObj += movement*cost[iSequence];
2209          matrix_->add(this,outputArray,iSequence,movement);
2210        } else if (value<tolerance) {
2211          // at correct bound but may swap
2212          FakeBound bound = getFakeBound(iSequence);
2213          if (bound==ClpSimplexDual::lowerFake) {
2214            movement = upper[iSequence]-lower[iSequence];
2215            assert (fabs(movement)<1.0e30);
2216            setStatus(iSequence,atUpperBound);
2217            solution[iSequence] = upper[iSequence];
2218            changeObj += movement*cost[iSequence];
2219            numberFake_--;
2220            setFakeBound(iSequence,noFake);
2221          }
2222        }
2223      } else if (status==atUpperBound) {
2224        double movement=0.0;
2225
2226        if (value>tolerance) {
2227          // to lower bound (if swap)
2228          // put back alpha
2229          which[numberInfeasibilities++]=iSequence;
2230          movement = lower[iSequence]-upper[iSequence];
2231          assert (fabs(movement)<1.0e30);
2232          changeObj += movement*cost[iSequence];
2233          matrix_->add(this,outputArray,iSequence,movement);
2234        } else if (value>-tolerance) {
2235          // at correct bound but may swap
2236          FakeBound bound = getFakeBound(iSequence);
2237          if (bound==ClpSimplexDual::upperFake) {
2238            movement = lower[iSequence]-upper[iSequence];
2239            assert (fabs(movement)<1.0e30);
2240            setStatus(iSequence,atLowerBound);
2241            solution[iSequence] = lower[iSequence];
2242            changeObj += movement*cost[iSequence];
2243            numberFake_--;
2244            setFakeBound(iSequence,noFake);
2245          }
2246        }
2247      }
2248    }
2249  }
2250#ifdef CLP_DEBUG
2251  if (fullRecompute&&numberFake_&&(handler_->logLevel()&16)!=0) 
2252    printf("%d fake after full update\n",numberFake_);
2253#endif
2254  // set number of infeasibilities
2255  columnArray->setNumElements(numberInfeasibilities);
2256  numberInfeasibilities += numberRowInfeasibilities;
2257  if (fullRecompute) {
2258    // do actual flips
2259    flipBounds(rowArray,columnArray,0.0);
2260  }
2261  objectiveChange += changeObj;
2262  return numberInfeasibilities;
2263}
2264void
2265ClpSimplexDual::updateDualsInValuesPass(CoinIndexedVector * rowArray,
2266                                  CoinIndexedVector * columnArray,
2267                                        double theta)
2268{
2269 
2270  // use a tighter tolerance except for all being okay
2271  double tolerance = dualTolerance_;
2272 
2273  // Coding is very similar but we can save a bit by splitting
2274  // Do rows
2275  {
2276    int i;
2277    double * reducedCost = djRegion(0);
2278    double * work;
2279    int number;
2280    int * which;
2281    work = rowArray->denseVector();
2282    number = rowArray->getNumElements();
2283    which = rowArray->getIndices();
2284    for (i=0;i<number;i++) {
2285      int iSequence = which[i];
2286      double alphaI = work[i];
2287      double value = reducedCost[iSequence]-theta*alphaI;
2288      work[i]=0.0;
2289      reducedCost[iSequence]=value;
2290     
2291      Status status = getStatus(iSequence+numberColumns_);
2292      // more likely to be at upper bound ?
2293      if (status==atUpperBound) {
2294
2295        if (value>tolerance) 
2296          reducedCost[iSequence]=0.0;
2297      } else if (status==atLowerBound) {
2298
2299        if (value<-tolerance) {
2300          reducedCost[iSequence]=0.0;
2301        }
2302      }
2303    }
2304  }
2305  rowArray->setNumElements(0);
2306
2307  // Do columns
2308  {
2309    int i;
2310    double * reducedCost = djRegion(1);
2311    double * work;
2312    int number;
2313    int * which;
2314    work = columnArray->denseVector();
2315    number = columnArray->getNumElements();
2316    which = columnArray->getIndices();
2317   
2318    for (i=0;i<number;i++) {
2319      int iSequence = which[i];
2320      double alphaI = work[i];
2321      double value = reducedCost[iSequence]-theta*alphaI;
2322      work[i]=0.0;
2323      reducedCost[iSequence]=value;
2324     
2325      Status status = getStatus(iSequence);
2326      if (status==atLowerBound) {
2327        if (value<-tolerance) 
2328          reducedCost[iSequence]=0.0;
2329      } else if (status==atUpperBound) {
2330        if (value>tolerance) 
2331          reducedCost[iSequence]=0.0;
2332      }
2333    }
2334  }
2335  columnArray->setNumElements(0);
2336}
2337/*
2338   Chooses dual pivot row
2339   Would be faster with separate region to scan
2340   and will have this (with square of infeasibility) when steepest
2341   For easy problems we can just choose one of the first rows we look at
2342*/
2343void
2344ClpSimplexDual::dualRow(int alreadyChosen)
2345{
2346  // get pivot row using whichever method it is
2347  int chosenRow=-1;
2348#ifdef FORCE_FOLLOW
2349  bool forceThis=false;
2350  if (!fpFollow&&strlen(forceFile)) {
2351    fpFollow = fopen(forceFile,"r");
2352    assert (fpFollow);
2353  }
2354  if (fpFollow) {
2355    if (numberIterations_<=force_iteration) {
2356      // read to next Clp0102
2357      char temp[300];
2358      while (fgets(temp,250,fpFollow)) {
2359        if (strncmp(temp,"Clp0102",7))
2360          continue;
2361        char cin,cout;
2362        sscanf(temp+9,"%d%*f%*s%*c%c%d%*s%*c%c%d",
2363               &force_iteration,&cin,&force_in,&cout,&force_out);
2364        if (cin=='R')
2365          force_in += numberColumns_;
2366        if (cout=='R')
2367          force_out += numberColumns_;
2368        forceThis=true;
2369        assert (numberIterations_==force_iteration-1);
2370        printf("Iteration %d will force %d out and %d in\n",
2371               force_iteration,force_out,force_in);
2372        alreadyChosen=force_out;
2373        break;
2374      }
2375    } else {
2376      // use old
2377      forceThis=true;
2378    }
2379    if (!forceThis) {
2380      fclose(fpFollow);
2381      fpFollow=NULL;
2382      forceFile="";
2383    }
2384  }
2385#endif
2386  double freeAlpha=0.0;
2387  if (alreadyChosen<0) {
2388    // first see if any free variables and put them in basis
2389    int nextFree = nextSuperBasic();
2390    //nextFree=-1; //off
2391    if (nextFree>=0) {
2392      // unpack vector and find a good pivot
2393      unpack(rowArray_[1],nextFree);
2394      factorization_->updateColumn(rowArray_[2],rowArray_[1]);
2395     
2396      double * work=rowArray_[1]->denseVector();
2397      int number=rowArray_[1]->getNumElements();
2398      int * which=rowArray_[1]->getIndices();
2399      double bestFeasibleAlpha=0.0;
2400      int bestFeasibleRow=-1;
2401      double bestInfeasibleAlpha=0.0;
2402      int bestInfeasibleRow=-1;
2403      int i;
2404     
2405      for (i=0;i<number;i++) {
2406        int iRow = which[i];
2407        double alpha = fabs(work[iRow]);
2408        if (alpha>1.0e-3) {
2409          int iSequence=pivotVariable_[iRow];
2410          double value = solution_[iSequence];
2411          double lower = lower_[iSequence];
2412          double upper = upper_[iSequence];
2413          double infeasibility=0.0;
2414          if (value>upper)
2415            infeasibility = value-upper;
2416          else if (value<lower)
2417            infeasibility = lower-value;
2418          if (infeasibility*alpha>bestInfeasibleAlpha&&alpha>1.0e-1) {
2419            if (!flagged(iSequence)) {
2420              bestInfeasibleAlpha = infeasibility*alpha;
2421              bestInfeasibleRow=iRow;
2422            }
2423          }
2424          if (alpha>bestFeasibleAlpha&&(lower>-1.0e20||upper<1.0e20)) {
2425            bestFeasibleAlpha = alpha;
2426            bestFeasibleRow=iRow;
2427          }
2428        }
2429      }
2430      if (bestInfeasibleRow>=0) 
2431        chosenRow = bestInfeasibleRow;
2432      else if (bestFeasibleAlpha>1.0e-2)
2433        chosenRow = bestFeasibleRow;
2434      if (chosenRow>=0) {
2435        pivotRow_=chosenRow;
2436        freeAlpha=work[chosenRow];
2437      }
2438      rowArray_[1]->clear();
2439    } 
2440  } else {
2441    // in values pass
2442    chosenRow=alreadyChosen;
2443#ifdef FORCE_FOLLOW
2444    if(forceThis) {
2445      alreadyChosen=-1;
2446      chosenRow=-1;
2447      for (int i=0;i<numberRows_;i++) {
2448        if (pivotVariable_[i]==force_out) {
2449          chosenRow=i;
2450          break;
2451        }
2452      }
2453      assert (chosenRow>=0);
2454    }
2455#endif
2456    pivotRow_=chosenRow;
2457  }
2458  if (chosenRow<0) 
2459    pivotRow_=dualRowPivot_->pivotRow();
2460
2461  if (pivotRow_>=0) {
2462    sequenceOut_ = pivotVariable_[pivotRow_];
2463    valueOut_ = solution_[sequenceOut_];
2464    lowerOut_ = lower_[sequenceOut_];
2465    upperOut_ = upper_[sequenceOut_];
2466    if (alreadyChosen<0) {
2467      // if we have problems we could try other way and hope we get a
2468      // zero pivot?
2469      if (valueOut_>upperOut_) {
2470        directionOut_ = -1;
2471        dualOut_ = valueOut_ - upperOut_;
2472      } else if (valueOut_<lowerOut_) {
2473        directionOut_ = 1;
2474        dualOut_ = lowerOut_ - valueOut_;
2475      } else {
2476#if 1
2477        // odd (could be free) - it's feasible - go to nearest
2478        if (valueOut_-lowerOut_<upperOut_-valueOut_) {
2479          directionOut_ = 1;
2480          dualOut_ = lowerOut_ - valueOut_;
2481        } else {
2482          directionOut_ = -1;
2483          dualOut_ = valueOut_ - upperOut_;
2484        }
2485#else
2486        // odd (could be free) - it's feasible - improve obj
2487        printf("direction from alpha of %g is %d\n",
2488               freeAlpha,freeAlpha>0.0 ? 1 : -1);
2489        if (valueOut_-lowerOut_>1.0e20)
2490          freeAlpha=1.0;
2491        else if(upperOut_-valueOut_>1.0e20) 
2492          freeAlpha=-1.0;
2493        //if (valueOut_-lowerOut_<upperOut_-valueOut_) {
2494        if (freeAlpha<0.0) {
2495          directionOut_ = 1;
2496          dualOut_ = lowerOut_ - valueOut_;
2497        } else {
2498          directionOut_ = -1;
2499          dualOut_ = valueOut_ - upperOut_;
2500        }
2501        printf("direction taken %d - bounds %g %g %g\n",
2502               directionOut_,lowerOut_,valueOut_,upperOut_);
2503#endif
2504      }
2505#ifdef CLP_DEBUG
2506      assert(dualOut_>=0.0);
2507#endif
2508    } else {
2509      // in values pass so just use sign of dj
2510      // We don't want to go through any barriers so set dualOut low
2511      // free variables will never be here
2512      dualOut_ = 1.0e-6;
2513      if (dj_[sequenceOut_]>0.0) {
2514        // this will give a -1 in pivot row (as slacks are -1.0)
2515        directionOut_ = 1;
2516      } else {
2517        directionOut_ = -1;
2518      }
2519    }
2520  }
2521  return ;
2522}
2523// Checks if any fake bounds active - if so returns number and modifies
2524// dualBound_ and everything.
2525// Free variables will be left as free
2526// Returns number of bounds changed if >=0
2527// Returns -1 if not initialize and no effect
2528// Fills in changeVector which can be used to see if unbounded
2529// and cost of change vector
2530int
2531ClpSimplexDual::changeBounds(int initialize,
2532                                 CoinIndexedVector * outputArray,
2533                                 double & changeCost)
2534{ 
2535  numberFake_ = 0;
2536  if (!initialize) {
2537    int numberInfeasibilities;
2538    double newBound;
2539    newBound = 5.0*dualBound_;
2540    numberInfeasibilities=0;
2541    changeCost=0.0;
2542    // put back original bounds and then check
2543    createRim1(false);
2544    int iSequence;
2545    // bounds will get bigger - just look at ones at bounds
2546    for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
2547      double lowerValue=lower_[iSequence];
2548      double upperValue=upper_[iSequence];
2549      double value=solution_[iSequence];
2550      setFakeBound(iSequence,ClpSimplexDual::noFake);
2551      switch(getStatus(iSequence)) {
2552
2553      case basic:
2554      case ClpSimplex::isFixed:
2555        break;
2556      case isFree:
2557      case superBasic:
2558        break;
2559      case atUpperBound:
2560        if (fabs(value-upperValue)>primalTolerance_) 
2561          numberInfeasibilities++;
2562        break;
2563      case atLowerBound:
2564        if (fabs(value-lowerValue)>primalTolerance_) 
2565          numberInfeasibilities++;
2566        break;
2567      }
2568    }
2569    // If dual infeasible then carry on
2570    if (numberInfeasibilities) {
2571      handler_->message(CLP_DUAL_CHECKB,messages_)
2572        <<newBound
2573        <<CoinMessageEol;
2574      int iSequence;
2575      for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
2576        double lowerValue=lower_[iSequence];
2577        double upperValue=upper_[iSequence];
2578        double newLowerValue;
2579        double newUpperValue;
2580        Status status = getStatus(iSequence);
2581        if (status==atUpperBound||
2582            status==atLowerBound) {
2583          double value = solution_[iSequence];
2584          if (value-lowerValue<=upperValue-value) {
2585            newLowerValue = CoinMax(lowerValue,value-0.666667*newBound);
2586            newUpperValue = CoinMin(upperValue,newLowerValue+newBound);
2587          } else {
2588            newUpperValue = CoinMin(upperValue,value+0.666667*newBound);
2589            newLowerValue = CoinMax(lowerValue,newUpperValue-newBound);
2590          }
2591          lower_[iSequence]=newLowerValue;
2592          upper_[iSequence]=newUpperValue;
2593          if (newLowerValue > lowerValue) {
2594            if (newUpperValue < upperValue) {
2595              setFakeBound(iSequence,ClpSimplexDual::bothFake);
2596              numberFake_++;
2597            } else {
2598              setFakeBound(iSequence,ClpSimplexDual::lowerFake);
2599              numberFake_++;
2600            }
2601          } else {
2602            if (newUpperValue < upperValue) {
2603              setFakeBound(iSequence,ClpSimplexDual::upperFake);
2604              numberFake_++;
2605            }
2606          }
2607          if (status==atUpperBound)
2608            solution_[iSequence] = newUpperValue;
2609          else
2610            solution_[iSequence] = newLowerValue;
2611          double movement = solution_[iSequence] - value;
2612          if (movement&&outputArray) {
2613            if (iSequence>=numberColumns_) {
2614              outputArray->quickAdd(iSequence,-movement);
2615              changeCost += movement*cost_[iSequence];
2616            } else {
2617              matrix_->add(this,outputArray,iSequence,movement);
2618              changeCost += movement*cost_[iSequence];
2619            }
2620          }
2621        }
2622      }
2623      dualBound_ = newBound;
2624    } else {
2625      numberInfeasibilities=-1;
2626    }
2627    return numberInfeasibilities;
2628  } else if (initialize==1) {
2629    int iSequence;
2630     
2631    for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
2632      Status status = getStatus(iSequence);
2633      if (status==atUpperBound||
2634          status==atLowerBound) {
2635        double lowerValue=lower_[iSequence];
2636        double upperValue=upper_[iSequence];
2637        double value = solution_[iSequence];
2638        if (lowerValue>-largeValue_||upperValue<largeValue_) {
2639          if (lowerValue-value>-0.5*dualBound_||
2640              upperValue-value<0.5*dualBound_) {
2641            if (fabs(lowerValue-value)<=fabs(upperValue-value)) {
2642              if (upperValue > lowerValue + dualBound_) {
2643                upper_[iSequence]=lowerValue+dualBound_;
2644                setFakeBound(iSequence,ClpSimplexDual::upperFake);
2645                numberFake_++;
2646              }
2647            } else {
2648              if (lowerValue < upperValue - dualBound_) {
2649                lower_[iSequence]=upperValue-dualBound_;
2650                setFakeBound(iSequence,ClpSimplexDual::lowerFake);
2651                numberFake_++;
2652              }
2653            }
2654          } else {
2655            lower_[iSequence]=-0.5*dualBound_;
2656            upper_[iSequence]= 0.5*dualBound_;
2657            setFakeBound(iSequence,ClpSimplexDual::bothFake);
2658            numberFake_++;
2659          }
2660          if (status==atUpperBound)
2661            solution_[iSequence]=upper_[iSequence];
2662          else
2663            solution_[iSequence]=lower_[iSequence];
2664        } else {
2665          // set non basic free variables to fake bounds
2666          // I don't think we should ever get here
2667          CoinAssert(!("should not be here"));
2668          lower_[iSequence]=-0.5*dualBound_;
2669          upper_[iSequence]= 0.5*dualBound_;
2670          setFakeBound(iSequence,ClpSimplexDual::bothFake);
2671          numberFake_++;
2672          setStatus(iSequence,atUpperBound);
2673          solution_[iSequence]=0.5*dualBound_;
2674        }
2675      }
2676    }
2677
2678    return 1;
2679  } else {
2680    // just reset changed ones
2681    if (columnScale_) {
2682      int iSequence;
2683      for (iSequence=0;iSequence<numberColumns_;iSequence++) {
2684        FakeBound fakeStatus = getFakeBound(iSequence);
2685        if (fakeStatus!=noFake) {
2686          if (((int) fakeStatus&1)!=0) {
2687            // lower
2688            double value = columnLower_[iSequence];
2689            if (value>-1.0e30) {
2690              double multiplier = rhsScale_/columnScale_[iSequence];
2691              value *= multiplier;
2692            }
2693            columnLowerWork_[iSequence]=value;
2694          }
2695          if (((int) fakeStatus&2)!=0) {
2696            // upper
2697            double value = columnUpper_[iSequence];
2698            if (value<1.0e30) {
2699              double multiplier = rhsScale_/columnScale_[iSequence];
2700              value *= multiplier;
2701            }
2702            columnUpperWork_[iSequence]=value;
2703          }
2704        }
2705      }
2706      for (iSequence=0;iSequence<numberRows_;iSequence++) {
2707        FakeBound fakeStatus = getFakeBound(iSequence+numberColumns_);
2708        if (fakeStatus!=noFake) {
2709          if (((int) fakeStatus&1)!=0) {
2710            // lower
2711            double value = rowLower_[iSequence];
2712            if (value>-1.0e30) {
2713              double multiplier = rhsScale_*rowScale_[iSequence];
2714              value *= multiplier;
2715            }
2716            rowLowerWork_[iSequence]=value;
2717          }
2718          if (((int) fakeStatus&2)!=0) {
2719            // upper
2720            double value = rowUpper_[iSequence];
2721            if (value<1.0e30) {
2722              double multiplier = rhsScale_*rowScale_[iSequence];
2723              value *= multiplier;
2724            }
2725            rowUpperWork_[iSequence]=value;
2726          }
2727        }
2728      }
2729    } else {
2730      int iSequence;
2731      for (iSequence=0;iSequence<numberColumns_;iSequence++) {
2732        FakeBound fakeStatus = getFakeBound(iSequence);
2733        if (((int) fakeStatus&1)!=0) {
2734          // lower
2735          columnLowerWork_[iSequence]=columnLower_[iSequence];
2736        }
2737        if (((int) fakeStatus&2)!=0) {
2738          // upper
2739          columnUpperWork_[iSequence]=columnUpper_[iSequence];
2740        }
2741      }
2742      for (iSequence=0;iSequence<numberRows_;iSequence++) {
2743        FakeBound fakeStatus = getFakeBound(iSequence+numberColumns_);
2744        if (((int) fakeStatus&1)!=0) {
2745          // lower
2746          rowLowerWork_[iSequence]=rowLower_[iSequence];
2747        }
2748        if (((int) fakeStatus&2)!=0) {
2749          // upper
2750          rowUpperWork_[iSequence]=rowUpper_[iSequence];
2751        }
2752      }
2753    }
2754    return 0;
2755  }
2756}
2757int
2758ClpSimplexDual::dualColumn0(const CoinIndexedVector * rowArray,
2759                           const CoinIndexedVector * columnArray,
2760                           CoinIndexedVector * spareArray,
2761                           double acceptablePivot,
2762                           double & upperReturn, double &bestReturn,double & badFree)
2763{
2764  // do first pass to get possibles
2765  double * spare = spareArray->denseVector();
2766  int * index = spareArray->getIndices();
2767  const double * work;
2768  int number;
2769  const int * which;
2770  const double * reducedCost;
2771  // We can also see if infeasible or pivoting on free
2772  double tentativeTheta = 1.0e25;
2773  double upperTheta = 1.0e31;
2774  double freePivot = acceptablePivot;
2775  double bestPossible=0.0;
2776  int numberRemaining=0;
2777  int i;
2778  badFree=0.0;
2779  for (int iSection=0;iSection<2;iSection++) {
2780
2781    int addSequence;
2782
2783    if (!iSection) {
2784      work = rowArray->denseVector();
2785      number = rowArray->getNumElements();
2786      which = rowArray->getIndices();
2787      reducedCost = rowReducedCost_;
2788      addSequence = numberColumns_;
2789    } else {
2790      work = columnArray->denseVector();
2791      number = columnArray->getNumElements();
2792      which = columnArray->getIndices();
2793      reducedCost = reducedCostWork_;
2794      addSequence = 0;
2795    }
2796   
2797    for (i=0;i<number;i++) {
2798      int iSequence = which[i];
2799      double alpha;
2800      double oldValue;
2801      double value;
2802      bool keep;
2803
2804      switch(getStatus(iSequence+addSequence)) {
2805         
2806      case basic:
2807      case ClpSimplex::isFixed:
2808        break;
2809      case isFree:
2810      case superBasic:
2811        alpha = work[i];
2812        bestPossible = CoinMax(bestPossible,fabs(alpha));
2813        oldValue = reducedCost[iSequence];
2814        // If free has to be very large - should come in via dualRow
2815        //if (getStatus(iSequence+addSequence)==isFree&&fabs(alpha)<1.0e-3)
2816        //break;
2817        if (oldValue>dualTolerance_) {
2818          keep = true;
2819        } else if (oldValue<-dualTolerance_) {
2820          keep = true;
2821        } else {
2822          if (fabs(alpha)>CoinMax(10.0*acceptablePivot,1.0e-5)) {
2823            keep = true;
2824          } else {
2825            keep=false;
2826            badFree=CoinMax(badFree,fabs(alpha));
2827          }
2828        }
2829        if (keep) {
2830          // free - choose largest
2831          if (fabs(alpha)>freePivot) {
2832            freePivot=fabs(alpha);
2833            sequenceIn_ = iSequence + addSequence;
2834            theta_=oldValue/alpha;
2835            alpha_=alpha;
2836          }
2837        }
2838        break;
2839      case atUpperBound:
2840        alpha = work[i];
2841        oldValue = reducedCost[iSequence];
2842        value = oldValue-tentativeTheta*alpha;
2843        //assert (oldValue<=dualTolerance_*1.0001);
2844        if (value>dualTolerance_) {
2845          bestPossible = CoinMax(bestPossible,-alpha);
2846          value = oldValue-upperTheta*alpha;
2847          if (value>dualTolerance_ && -alpha>=acceptablePivot) 
2848            upperTheta = (oldValue-dualTolerance_)/alpha;
2849          // add to list
2850          spare[numberRemaining]=alpha;
2851          index[numberRemaining++]=iSequence+addSequence;
2852        }
2853        break;
2854      case atLowerBound:
2855        alpha = work[i];
2856        oldValue = reducedCost[iSequence];
2857        value = oldValue-tentativeTheta*alpha;
2858        //assert (oldValue>=-dualTolerance_*1.0001);
2859        if (value<-dualTolerance_) {
2860          bestPossible = CoinMax(bestPossible,alpha);
2861          value = oldValue-upperTheta*alpha;
2862          if (value<-dualTolerance_ && alpha>=acceptablePivot) 
2863            upperTheta = (oldValue+dualTolerance_)/alpha;
2864          // add to list
2865          spare[numberRemaining]=alpha;
2866          index[numberRemaining++]=iSequence+addSequence;
2867        }
2868        break;
2869      }
2870    }
2871  }
2872  upperReturn = upperTheta;
2873  bestReturn = bestPossible;
2874  return numberRemaining;
2875}
2876/*
2877   Row array has row part of pivot row (as duals so sign may be switched)
2878   Column array has column part.
2879   This chooses pivot column.
2880   Spare array will be needed when we start getting clever.
2881   We will check for basic so spare array will never overflow.
2882   If necessary will modify costs
2883*/
2884double
2885ClpSimplexDual::dualColumn(CoinIndexedVector * rowArray,
2886                           CoinIndexedVector * columnArray,
2887                           CoinIndexedVector * spareArray,
2888                           CoinIndexedVector * spareArray2,
2889                           double acceptablePivot,
2890                           CoinBigIndex * dubiousWeights)
2891{
2892  int numberPossiblySwapped=0;
2893  int numberRemaining=0;
2894 
2895  double totalThru=0.0; // for when variables flip
2896  //double saveAcceptable=acceptablePivot;
2897  //acceptablePivot=1.0e-9;
2898
2899  double bestEverPivot=acceptablePivot;
2900  int lastSequence = -1;
2901  double lastPivot=0.0;
2902  double upperTheta;
2903  double newTolerance = dualTolerance_;
2904  //newTolerance = dualTolerance_+1.0e-6*dblParam_[ClpDualTolerance];
2905  // will we need to increase tolerance
2906  bool thisIncrease=false;
2907  // If we think we need to modify costs (not if something from broad sweep)
2908  bool modifyCosts=false;
2909  // Increase in objective due to swapping bounds (may be negative)
2910  double increaseInObjective=0.0;
2911
2912  // use spareArrays to put ones looked at in
2913  // we are going to flip flop between
2914  int iFlip = 0;
2915  // Possible list of pivots
2916  int interesting[2];
2917  // where possible swapped ones are
2918  int swapped[2];
2919  // for zeroing out arrays after
2920  int marker[2][2];
2921  // pivot elements
2922  double * array[2], * spare, * spare2;
2923  // indices
2924  int * indices[2], * index, * index2;
2925  spareArray2->clear();
2926  array[0] = spareArray->denseVector();
2927  indices[0] = spareArray->getIndices();
2928  spare = array[0];
2929  index = indices[0];
2930  array[1] = spareArray2->denseVector();
2931  indices[1] = spareArray2->getIndices();
2932  int i;
2933
2934  // initialize lists
2935  for (i=0;i<2;i++) {
2936    interesting[i]=0;
2937    swapped[i]=numberColumns_;
2938    marker[i][0]=0;
2939    marker[i][1]=numberColumns_;
2940  }
2941  /*
2942    First we get a list of possible pivots.  We can also see if the
2943    problem looks infeasible or whether we want to pivot in free variable.
2944    This may make objective go backwards but can only happen a finite
2945    number of times and I do want free variables basic.
2946
2947    Then we flip back and forth.  At the start of each iteration
2948    interesting[iFlip] should have possible candidates and swapped[iFlip]
2949    will have pivots if we decide to take a previous pivot.
2950    At end of each iteration interesting[1-iFlip] should have
2951    candidates if we go through this theta and swapped[1-iFlip]
2952    pivots if we don't go through.
2953
2954    At first we increase theta and see what happens.  We start
2955    theta at a reasonable guess.  If in right area then we do bit by bit.
2956
2957   */
2958
2959  // do first pass to get possibles
2960  upperTheta = 1.0e31;
2961  double bestPossible=0.0;
2962  double badFree=0.0;
2963  alpha_=0.0;
2964  if (spareIntArray_[0]!=-1) {
2965    numberRemaining = dualColumn0(rowArray,columnArray,spareArray,
2966                                  acceptablePivot,upperTheta,bestPossible,badFree);
2967  } else {
2968    // already done
2969    numberRemaining = spareArray->getNumElements();
2970    spareArray->setNumElements(0);
2971    upperTheta = spareDoubleArray_[0];
2972    bestPossible = spareDoubleArray_[1];
2973    theta_ = spareDoubleArray_[2];
2974    alpha_ = spareDoubleArray_[3];
2975    sequenceIn_ = spareIntArray_[1];
2976  }
2977  // switch off
2978  spareIntArray_[0]=0;
2979  // We can also see if infeasible or pivoting on free
2980  double tentativeTheta = 1.0e25;
2981  interesting[0]=numberRemaining;
2982  marker[0][0] = numberRemaining;
2983
2984  if (!numberRemaining&&sequenceIn_<0)
2985    return 0.0; // Looks infeasible
2986
2987  // If sum of bad small pivots too much
2988#define MORE_CAREFUL
2989#ifdef MORE_CAREFUL
2990  bool badSumPivots=false;
2991#endif
2992  if (sequenceIn_>=0) {
2993    // free variable - always choose
2994  } else {
2995
2996    theta_=1.0e50;
2997    // now flip flop between spare arrays until reasonable theta
2998    tentativeTheta = CoinMax(10.0*upperTheta,1.0e-7);
2999   
3000    // loops increasing tentative theta until can't go through
3001   
3002    while (tentativeTheta < 1.0e22) {
3003      double thruThis = 0.0;
3004     
3005      double bestPivot=acceptablePivot;
3006      int bestSequence=-1;
3007     
3008      numberPossiblySwapped = numberColumns_;
3009      numberRemaining = 0;
3010
3011      upperTheta = 1.0e50;
3012
3013      spare = array[iFlip];
3014      index = indices[iFlip];
3015      spare2 = array[1-iFlip];
3016      index2 = indices[1-iFlip];
3017
3018      // try 3 different ways
3019      // 1 bias increase by ones with slightly wrong djs
3020      // 2 bias by all
3021      // 3 bias by all - tolerance
3022#define TRYBIAS 3
3023
3024
3025      double increaseInThis=0.0; //objective increase in this loop
3026     
3027      for (i=0;i<interesting[iFlip];i++) {
3028        int iSequence = index[i];
3029        double alpha = spare[i];
3030        double oldValue = dj_[iSequence];
3031        double value = oldValue-tentativeTheta*alpha;
3032
3033        if (alpha < 0.0) {
3034          //at upper bound
3035          if (value>newTolerance) {
3036            double range = upper_[iSequence] - lower_[iSequence];
3037            thruThis -= range*alpha;
3038#if TRYBIAS==1
3039            if (oldValue>0.0)
3040              increaseInThis -= oldValue*range;
3041#elif TRYBIAS==2
3042            increaseInThis -= oldValue*range;
3043#else
3044            increaseInThis -= (oldValue+dualTolerance_)*range;
3045#endif
3046            // goes on swapped list (also means candidates if too many)
3047            spare2[--numberPossiblySwapped]=alpha;
3048            index2[numberPossiblySwapped]=iSequence;
3049            if (fabs(alpha)>bestPivot) {
3050              bestPivot=fabs(alpha);
3051              bestSequence=numberPossiblySwapped;
3052            }
3053          } else {
3054            value = oldValue-upperTheta*alpha;
3055            if (value>newTolerance && -alpha>=acceptablePivot) 
3056              upperTheta = (oldValue-newTolerance)/alpha;
3057            spare2[numberRemaining]=alpha;
3058            index2[numberRemaining++]=iSequence;
3059          }
3060        } else {
3061          // at lower bound
3062          if (value<-newTolerance) {
3063            double range = upper_[iSequence] - lower_[iSequence];
3064            thruThis += range*alpha;
3065            //?? is this correct - and should we look at good ones
3066#if TRYBIAS==1
3067            if (oldValue<0.0)
3068              increaseInThis += oldValue*range;
3069#elif TRYBIAS==2
3070            increaseInThis += oldValue*range;
3071#else
3072            increaseInThis += (oldValue-dualTolerance_)*range;
3073#endif
3074            // goes on swapped list (also means candidates if too many)
3075            spare2[--numberPossiblySwapped]=alpha;
3076            index2[numberPossiblySwapped]=iSequence;
3077            if (fabs(alpha)>bestPivot) {
3078              bestPivot=fabs(alpha);
3079              bestSequence=numberPossiblySwapped;
3080            }
3081          } else {
3082            value = oldValue-upperTheta*alpha;
3083            if (value<-newTolerance && alpha>=acceptablePivot) 
3084              upperTheta = (oldValue+newTolerance)/alpha;
3085            spare2[numberRemaining]=alpha;
3086            index2[numberRemaining++]=iSequence;
3087          }
3088        }
3089      }
3090      swapped[1-iFlip]=numberPossiblySwapped;
3091      interesting[1-iFlip]=numberRemaining;
3092      marker[1-iFlip][0]= CoinMax(marker[1-iFlip][0],numberRemaining);
3093      marker[1-iFlip][1]= CoinMin(marker[1-iFlip][1],numberPossiblySwapped);
3094     
3095      if (totalThru+thruThis>=fabs(dualOut_)||
3096          increaseInObjective+increaseInThis<0.0) {
3097        // We should be pivoting in this batch
3098        // so compress down to this lot
3099        numberRemaining=0;
3100        for (i=numberColumns_-1;i>=swapped[1-iFlip];i--) {
3101          spare[numberRemaining]=spare2[i];
3102          index[numberRemaining++]=index2[i];
3103        }
3104        interesting[iFlip]=numberRemaining;
3105        int iTry;
3106#define MAXTRY 100
3107        // first get ratio with tolerance
3108        for (iTry=0;iTry<MAXTRY;iTry++) {
3109         
3110          upperTheta=1.0e50;
3111          numberPossiblySwapped = numberColumns_;
3112          numberRemaining = 0;
3113
3114          increaseInThis=0.0; //objective increase in this loop
3115
3116          thruThis=0.0;
3117
3118          spare = array[iFlip];
3119          index = indices[iFlip];
3120          spare2 = array[1-iFlip];
3121          index2 = indices[1-iFlip];
3122          for (i=0;i<interesting[iFlip];i++) {
3123            int iSequence=index[i];
3124            double alpha=spare[i];
3125            double oldValue = dj_[iSequence];
3126            double value = oldValue-upperTheta*alpha;
3127           
3128            if (alpha < 0.0) {
3129              //at upper bound
3130              if (value>newTolerance) {
3131                if (-alpha>=acceptablePivot) {
3132                  upperTheta = (oldValue-newTolerance)/alpha;
3133                }
3134              }
3135            } else {
3136              // at lower bound
3137              if (value<-newTolerance) {
3138                if (alpha>=acceptablePivot) {
3139                  upperTheta = (oldValue+newTolerance)/alpha;
3140                }
3141              }
3142            }
3143          }
3144          bestPivot=acceptablePivot;
3145          sequenceIn_=-1;
3146#ifdef DUBIOUS_WEIGHTS
3147          double bestWeight=COIN_DBL_MAX;
3148#endif
3149          double largestPivot=acceptablePivot;
3150          // now choose largest and sum all ones which will go through
3151          //printf("XX it %d number %d\n",numberIterations_,interesting[iFlip]);
3152  // Sum of bad small pivots
3153#ifdef MORE_CAREFUL
3154          double sumBadPivots=0.0;
3155          badSumPivots=false;
3156#endif
3157          // Make sure upperTheta will work (-O2 and above gives problems)
3158          upperTheta *= 1.0000000001;
3159          for (i=0;i<interesting[iFlip];i++) {
3160            int iSequence=index[i];
3161            double alpha=spare[i];
3162            double value = dj_[iSequence]-upperTheta*alpha;
3163            double badDj=0.0;
3164           
3165            bool addToSwapped=false;
3166           
3167            if (alpha < 0.0) {
3168              //at upper bound
3169              if (value>=0.0) { 
3170                addToSwapped=true;
3171#if TRYBIAS==1
3172                badDj = -CoinMax(dj_[iSequence],0.0);
3173#elif TRYBIAS==2
3174                badDj = -dj_[iSequence];
3175#else
3176                badDj = -dj_[iSequence]-dualTolerance_;
3177#endif
3178              }
3179            } else {
3180              // at lower bound
3181              if (value<=0.0) {
3182                addToSwapped=true;
3183#if TRYBIAS==1
3184                badDj = CoinMin(dj_[iSequence],0.0);
3185#elif TRYBIAS==2
3186                badDj = dj_[iSequence];
3187#else
3188                badDj = dj_[iSequence]-dualTolerance_;
3189#endif
3190              }
3191            }
3192            if (!addToSwapped) {
3193              // add to list of remaining
3194              spare2[numberRemaining]=alpha;
3195              index2[numberRemaining++]=iSequence;
3196            } else {
3197              // add to list of swapped
3198              spare2[--numberPossiblySwapped]=alpha;
3199              index2[numberPossiblySwapped]=iSequence;
3200              // select if largest pivot
3201              bool take=false;
3202              double absAlpha = fabs(alpha);
3203#ifdef DUBIOUS_WEIGHTS
3204              // User could do anything to break ties here
3205              double weight;
3206              if (dubiousWeights)
3207                weight=dubiousWeights[iSequence];
3208              else
3209                weight=1.0;
3210              weight += randomNumberGenerator_.randomDouble()*1.0e-2;
3211              if (absAlpha>2.0*bestPivot) {
3212                take=true;
3213              } else if (absAlpha>largestPivot) {
3214                // could multiply absAlpha and weight
3215                if (weight*bestPivot<bestWeight*absAlpha)
3216                  take=true;
3217              }
3218#else
3219              if (absAlpha>bestPivot) 
3220                take=true;
3221#endif
3222#ifdef MORE_CAREFUL
3223              if (absAlpha<acceptablePivot&&upperTheta<1.0e20) {
3224                if (alpha < 0.0) {
3225                  //at upper bound
3226                  if (value>dualTolerance_) {
3227                    double gap=upper_[iSequence]-lower_[iSequence];
3228                    if (gap<1.0e20)
3229                      sumBadPivots += value*gap; 
3230                    else
3231                      sumBadPivots += 1.0e20;
3232                    //printf("bad %d alpha %g dj at upper %g\n",
3233                    //     iSequence,alpha,value);
3234                  }
3235                } else {
3236                  //at lower bound
3237                  if (value<-dualTolerance_) {
3238                    double gap=upper_[iSequence]-lower_[iSequence];
3239                    if (gap<1.0e20)
3240                      sumBadPivots -= value*gap; 
3241                    else
3242                      sumBadPivots += 1.0e20;
3243                    //printf("bad %d alpha %g dj at lower %g\n",
3244                    //     iSequence,alpha,value);
3245                  }
3246                }
3247              }
3248#endif
3249#ifdef FORCE_FOLLOW
3250              if (iSequence==force_in) {
3251                printf("taking %d - alpha %g best %g\n",force_in,absAlpha,largestPivot);
3252                take=true;
3253              }
3254#endif
3255              if (take) {
3256                sequenceIn_ = numberPossiblySwapped;
3257                bestPivot =  absAlpha;
3258                theta_ = dj_[iSequence]/alpha;
3259                largestPivot = CoinMax(largestPivot,0.5*bestPivot);
3260#ifdef DUBIOUS_WEIGHTS
3261                bestWeight = weight;
3262#endif
3263                //printf(" taken seq %d alpha %g weight %d\n",
3264                //   iSequence,absAlpha,dubiousWeights[iSequence]);
3265              } else {
3266                //printf(" not taken seq %d alpha %g weight %d\n",
3267                //   iSequence,absAlpha,dubiousWeights[iSequence]);
3268              }
3269              double range = upper_[iSequence] - lower_[iSequence];
3270              thruThis += range*fabs(alpha);
3271              increaseInThis += badDj*range;
3272            }
3273          }
3274#ifdef MORE_CAREFUL
3275          // If we have done pivots and things look bad set alpha_ 0.0 to force factorization
3276          if (sumBadPivots>1.0e4) {
3277            if (handler_->logLevel()>1)
3278            printf("maybe forcing re-factorization - sum %g  %d pivots\n",sumBadPivots,
3279                   factorization_->pivots());
3280            badSumPivots=true;
3281          }
3282#endif
3283          swapped[1-iFlip]=numberPossiblySwapped;
3284          interesting[1-iFlip]=numberRemaining;
3285          marker[1-iFlip][0]= CoinMax(marker[1-iFlip][0],numberRemaining);
3286          marker[1-iFlip][1]= CoinMin(marker[1-iFlip][1],numberPossiblySwapped);
3287          // If we stop now this will be increase in objective (I think)
3288          double increase = (fabs(dualOut_)-totalThru)*theta_;
3289          increase += increaseInObjective;
3290          if (theta_<0.0)
3291            thruThis += fabs(dualOut_); // force using this one
3292          if (increaseInObjective<0.0&&increase<0.0&&lastSequence>=0) {
3293            // back
3294            // We may need to be more careful - we could do by
3295            // switch so we always do fine grained?
3296            bestPivot=0.0;
3297          } else {
3298            // add in
3299            totalThru += thruThis;
3300            increaseInObjective += increaseInThis;
3301          }
3302          if (bestPivot<0.1*bestEverPivot&&
3303              bestEverPivot>1.0e-6&&
3304              (bestPivot<1.0e-3||totalThru*2.0>fabs(dualOut_))) {
3305            // back to previous one
3306            sequenceIn_=lastSequence;
3307            // swap regions
3308            iFlip = 1-iFlip;
3309            break;
3310          } else if (sequenceIn_==-1&&upperTheta>largeValue_) {
3311            if (lastPivot>acceptablePivot) {
3312              // back to previous one
3313              sequenceIn_=lastSequence;
3314              // swap regions
3315              iFlip = 1-iFlip;
3316            } else {
3317              // can only get here if all pivots too small
3318            }
3319            break;
3320          } else if (totalThru>=fabs(dualOut_)) {
3321            modifyCosts=true; // fine grain - we can modify costs
3322            break; // no point trying another loop
3323          } else {
3324            lastSequence=sequenceIn_;
3325            if (bestPivot>bestEverPivot)
3326              bestEverPivot=bestPivot;
3327            iFlip = 1 -iFlip;
3328            modifyCosts=true; // fine grain - we can modify costs
3329          }
3330        }
3331        if (iTry==MAXTRY)
3332          iFlip = 1-iFlip; // flip back
3333        break;
3334      } else {
3335        // skip this lot
3336        if (bestPivot>1.0e-3||bestPivot>bestEverPivot) {
3337          bestEverPivot=bestPivot;
3338          lastSequence=bestSequence;
3339        } else {
3340          // keep old swapped
3341          CoinMemcpyN(array[iFlip]+swapped[iFlip],
3342                 numberColumns_-swapped[iFlip],array[1-iFlip]+swapped[iFlip]);
3343          CoinMemcpyN(indices[iFlip]+swapped[iFlip],
3344                 numberColumns_-swapped[iFlip],indices[1-iFlip]+swapped[iFlip]);
3345          marker[1-iFlip][1] = CoinMin(marker[1-iFlip][1],swapped[iFlip]);
3346          swapped[1-iFlip]=swapped[iFlip];
3347        }
3348        increaseInObjective += increaseInThis;
3349        iFlip = 1 - iFlip; // swap regions
3350        tentativeTheta = 2.0*upperTheta;
3351        totalThru += thruThis;
3352      }
3353    }
3354   
3355    // can get here without sequenceIn_ set but with lastSequence
3356    if (sequenceIn_<0&&lastSequence>=0) {
3357      // back to previous one
3358      sequenceIn_=lastSequence;
3359      // swap regions
3360      iFlip = 1-iFlip;
3361    }
3362   
3363#define MINIMUMTHETA 1.0e-18
3364    // Movement should be minimum for anti-degeneracy - unless
3365    // fixed variable out
3366    double minimumTheta;
3367    if (upperOut_>lowerOut_)
3368      minimumTheta=MINIMUMTHETA;
3369    else
3370      minimumTheta=0.0;
3371    if (sequenceIn_>=0) {
3372      // at this stage sequenceIn_ is just pointer into index array
3373      // flip just so we can use iFlip
3374      iFlip = 1 -iFlip;
3375      spare = array[iFlip];
3376      index = indices[iFlip];
3377      double oldValue;
3378      alpha_ = spare[sequenceIn_];
3379      sequenceIn_ = indices[iFlip][sequenceIn_];
3380      oldValue = dj_[sequenceIn_];
3381      theta_ = CoinMax(oldValue/alpha_,0.0);
3382      if (theta_<minimumTheta&&fabs(alpha_)<1.0e5&&1) {
3383        // can't pivot to zero
3384#if 0
3385        if (oldValue-minimumTheta*alpha_>=-dualTolerance_) {
3386          theta_=minimumTheta;
3387        } else if (oldValue-minimumTheta*alpha_>=-newTolerance) {
3388          theta_=minimumTheta;
3389          thisIncrease=true;
3390        } else {
3391          theta_=CoinMax((oldValue+newTolerance)/alpha_,0.0);
3392          thisIncrease=true;
3393        }
3394#else
3395        theta_=minimumTheta;
3396#endif
3397      }
3398      // may need to adjust costs so all dual feasible AND pivoted is exactly 0
3399      //int costOffset = numberRows_+numberColumns_;
3400      if (modifyCosts) {
3401        int i;
3402        for (i=numberColumns_-1;i>=swapped[iFlip];i--) {
3403          int iSequence=index[i];
3404          double alpha=spare[i];
3405          double value = dj_[iSequence]-theta_*alpha;
3406           
3407          // can't be free here
3408         
3409          if (alpha < 0.0) {
3410            //at upper bound
3411            if (value>dualTolerance_) {
3412              thisIncrease=true;
3413#define MODIFYCOST 2
3414#if MODIFYCOST
3415              // modify cost to hit new tolerance
3416              double modification = alpha*theta_-dj_[iSequence]
3417                +newTolerance;
3418              if ((specialOptions_&(2048+4096+16384))!=0) {
3419                if ((specialOptions_&16384)!=0) {
3420                  if (fabs(modification)<1.0e-8)
3421                    modification=0.0;
3422                } else if ((specialOptions_&2048)!=0) {
3423                  if (fabs(modification)<1.0e-10)
3424                    modification=0.0;
3425                } else {
3426                  if (fabs(modification)<1.0e-12)
3427                    modification=0.0;
3428                }
3429              }
3430              dj_[iSequence] += modification;
3431              cost_[iSequence] +=  modification;
3432              if (modification)
3433                numberChanged_ ++; // Say changed costs
3434              //cost_[iSequence+costOffset] += modification; // save change
3435#endif
3436            }
3437          } else {
3438            // at lower bound
3439            if (-value>dualTolerance_) {
3440              thisIncrease=true;
3441#if MODIFYCOST
3442              // modify cost to hit new tolerance
3443              double modification = alpha*theta_-dj_[iSequence]
3444                -newTolerance;
3445              //modification = CoinMax(modification,-dualTolerance_);
3446              //assert (fabs(modification)<1.0e-7);
3447              if ((specialOptions_&(2048+4096))!=0) {
3448                if ((specialOptions_&2048)!=0) {
3449                  if (fabs(modification)<1.0e-10)
3450                    modification=0.0;
3451                } else {
3452                  if (fabs(modification)<1.0e-12)
3453                    modification=0.0;
3454                }
3455              }
3456              dj_[iSequence] += modification;
3457              cost_[iSequence] +=  modification;
3458              if (modification)
3459                numberChanged_ ++; // Say changed costs
3460              //cost_[iSequence+costOffset] += modification; // save change
3461#endif
3462            }
3463          }
3464        }
3465      }
3466    }
3467  }
3468
3469  if (sequenceIn_>=0) {
3470#ifdef MORE_CAREFUL
3471    // If we have done pivots and things look bad set alpha_ 0.0 to force factorization
3472    if (badSumPivots&&factorization_->pivots()) {
3473      if (handler_->logLevel()>1)
3474        printf("forcing re-factorization\n");
3475      alpha_=0.0;
3476    }
3477    if (fabs(theta_*badFree)>10.0*dualTolerance_&&factorization_->pivots()) {
3478      if (handler_->logLevel()>1)
3479        printf("forcing re-factorizationon free\n");
3480      alpha_=0.0;
3481    }
3482#endif
3483    lowerIn_ = lower_[sequenceIn_];
3484    upperIn_ = upper_[sequenceIn_];
3485    valueIn_ = solution_[sequenceIn_];
3486    dualIn_ = dj_[sequenceIn_];
3487
3488    if (numberTimesOptimal_) {
3489      // can we adjust cost back closer to original
3490      //*** add coding
3491    }
3492#if MODIFYCOST>1   
3493    // modify cost to hit zero exactly
3494    // so (dualIn_+modification)==theta_*alpha_
3495    double modification = theta_*alpha_-dualIn_;
3496    // But should not move objectivetoo much ??
3497#define DONT_MOVE_OBJECTIVE
3498#ifdef DONT_MOVE_OBJECTIVE
3499    double moveObjective = fabs(modification*solution_[sequenceIn_]);
3500    double smallMove = CoinMax(fabs(objectiveValue_),1.0e-3);
3501    if (moveObjective>smallMove) {
3502      if (handler_->logLevel()>1)
3503        printf("would move objective by %g - original mod %g sol value %g\n",moveObjective,
3504               modification,solution_[sequenceIn_]);
3505      modification *= smallMove/moveObjective;
3506    }
3507#endif
3508    if ((specialOptions_&(2048+4096))!=0) {
3509      if ((specialOptions_&16384)!=0) {
3510        // in fast dual
3511        if (fabs(modification)<1.0e-7)
3512          modification=0.0;
3513      } else if ((specialOptions_&2048)!=0) {
3514        if (fabs(modification)<1.0e-10)
3515          modification=0.0;
3516      } else {
3517        if (fabs(modification)<1.0e-12)
3518          modification=0.0;
3519      }
3520    }
3521    dualIn_ += modification;
3522    dj_[sequenceIn_]=dualIn_;
3523    cost_[sequenceIn_] += modification;
3524    if (modification)
3525      numberChanged_ ++; // Say changed costs
3526    //int costOffset = numberRows_+numberColumns_;
3527    //cost_[sequenceIn_+costOffset] += modification; // save change
3528    //assert (fabs(modification)<1.0e-6);
3529#ifdef CLP_DEBUG
3530    if ((handler_->logLevel()&32)&&fabs(modification)>1.0e-15)
3531      printf("exact %d new cost %g, change %g\n",sequenceIn_,
3532             cost_[sequenceIn_],modification);
3533#endif
3534#endif
3535   
3536    if (alpha_<0.0) {
3537      // as if from upper bound
3538      directionIn_=-1;
3539      upperIn_=valueIn_;
3540    } else {
3541      // as if from lower bound
3542      directionIn_=1;
3543      lowerIn_=valueIn_;
3544    }
3545  }
3546  //if (thisIncrease)
3547  //dualTolerance_+= 1.0e-6*dblParam_[ClpDualTolerance];
3548
3549  // clear arrays
3550
3551  for (i=0;i<2;i++) {
3552    CoinZeroN(array[i],marker[i][0]);
3553    CoinZeroN(array[i]+marker[i][1],numberColumns_-marker[i][1]);
3554  }
3555  return bestPossible;
3556}
3557#ifdef CLP_ALL_ONE_FILE
3558#undef MAXTRY
3559#endif
3560/* Checks if tentative optimal actually means unbounded
3561   Returns -3 if not, 2 if is unbounded */
3562int 
3563ClpSimplexDual::checkUnbounded(CoinIndexedVector * ray,
3564                                   CoinIndexedVector * spare,
3565                                   double changeCost)
3566{
3567  int status=2; // say unbounded
3568  factorization_->updateColumn(spare,ray);
3569  // get reduced cost
3570  int i;
3571  int number=ray->getNumElements();
3572  int * index = ray->getIndices();
3573  double * array = ray->denseVector();
3574  for (i=0;i<number;i++) {
3575    int iRow=index[i];
3576    int iPivot=pivotVariable_[iRow];
3577    changeCost -= cost(iPivot)*array[iRow];
3578  }
3579  double way;
3580  if (changeCost>0.0) {
3581    //try going down
3582    way=1.0;
3583  } else if (changeCost<0.0) {
3584    //try going up
3585    way=-1.0;
3586  } else {
3587#ifdef CLP_DEBUG
3588    printf("can't decide on up or down\n");
3589#endif
3590    way=0.0;
3591    status=-3;
3592  }
3593  double movement=1.0e10*way; // some largish number
3594  double zeroTolerance = 1.0e-14*dualBound_;
3595  for (i=0;i<number;i++) {
3596    int iRow=index[i];
3597    int iPivot=pivotVariable_[iRow];
3598    double arrayValue = array[iRow];
3599    if (fabs(arrayValue)<zeroTolerance)
3600      arrayValue=0.0;
3601    double newValue=solution(iPivot)+movement*arrayValue;
3602    if (newValue>upper(iPivot)+primalTolerance_||
3603        newValue<lower(iPivot)-primalTolerance_)
3604      status=-3; // not unbounded
3605  }
3606  if (status==2) {
3607    // create ray
3608    delete [] ray_;
3609    ray_ = new double [numberColumns_];
3610    CoinZeroN(ray_,numberColumns_);
3611    for (i=0;i<number;i++) {
3612      int iRow=index[i];
3613      int iPivot=pivotVariable_[iRow];
3614      double arrayValue = array[iRow];
3615      if (iPivot<numberColumns_&&fabs(arrayValue)>=zeroTolerance)
3616        ray_[iPivot] = way* array[iRow];
3617    }
3618  }
3619  ray->clear();
3620  return status;
3621}
3622//static int count_alpha=0;
3623/* Checks if finished.  Updates status */
3624void 
3625ClpSimplexDual::statusOfProblemInDual(int & lastCleaned,int type,
3626                                      double * givenDuals, ClpDataSave & saveData,
3627                                      int ifValuesPass)
3628{
3629#ifdef CLP_INVESTIGATE
3630  if (z_thinks>0&&z_thinks<2)
3631    z_thinks+=2;
3632#endif
3633  // If lots of iterations then adjust costs if large ones
3634  if (numberIterations_>4*(numberRows_+numberColumns_)&&objectiveScale_==1.0) {
3635    double largest=0.0;
3636    for (int i=0;i<numberRows_;i++) {
3637      int iColumn = pivotVariable_[i];
3638      largest = CoinMax(largest,fabs(cost_[iColumn]));
3639    }
3640    if (largest>1.0e6) {
3641      objectiveScale_ = 1.0e6/largest;
3642      for (int i=0;i<numberRows_+numberColumns_;i++)
3643        cost_[i] *= objectiveScale_;
3644    }
3645  }
3646  int numberPivots = factorization_->pivots();
3647  double realDualInfeasibilities=0.0;
3648  if (type==2) {
3649    if (alphaAccuracy_!=-1.0)
3650      alphaAccuracy_=-2.0;
3651    // trouble - restore solution
3652    CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
3653    CoinMemcpyN(savedSolution_+numberColumns_ ,
3654           numberRows_,rowActivityWork_);
3655    CoinMemcpyN(savedSolution_ ,
3656           numberColumns_,columnActivityWork_);
3657    // restore extra stuff
3658    int dummy;
3659    matrix_->generalExpanded(this,6,dummy);
3660    forceFactorization_=1; // a bit drastic but ..
3661    changeMade_++; // say something changed
3662    // get correct bounds on all variables
3663    resetFakeBounds();
3664  }
3665  int tentativeStatus = problemStatus_;
3666  double changeCost;
3667  bool unflagVariables = true;
3668  bool weightsSaved=false;
3669  int dontFactorizePivots = dontFactorizePivots_;
3670  if (type==3) {
3671    type=1;
3672    dontFactorizePivots=1;
3673  }
3674  if (alphaAccuracy_<0.0||!numberPivots||alphaAccuracy_>1.0e4||numberPivots>20) {
3675    if (problemStatus_>-3||numberPivots>dontFactorizePivots) {
3676      // factorize
3677      // later on we will need to recover from singularities
3678      // also we could skip if first time
3679      // save dual weights
3680      dualRowPivot_->saveWeights(this,1);
3681      weightsSaved=true;
3682      if (type) {
3683        // is factorization okay?
3684        if (internalFactorize(1)) {
3685          // no - restore previous basis
3686          unflagVariables = false;
3687          assert (type==1);
3688          changeMade_++; // say something changed
3689          // Keep any flagged variables
3690          int i;
3691          for (i=0;i<numberRows_+numberColumns_;i++) {
3692            if (flagged(i))
3693              saveStatus_[i] |= 64; //say flagged
3694          }
3695          CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
3696          CoinMemcpyN(savedSolution_+numberColumns_ ,
3697                      numberRows_,rowActivityWork_);
3698          CoinMemcpyN(savedSolution_ ,
3699                      numberColumns_,columnActivityWork_);
3700          // restore extra stuff
3701          int dummy;
3702          matrix_->generalExpanded(this,6,dummy);
3703          // get correct bounds on all variables
3704          resetFakeBounds();
3705          // need to reject something
3706          char x = isColumn(sequenceOut_) ? 'C' :'R';
3707          handler_->message(CLP_SIMPLEX_FLAG,messages_)
3708            <<x<<sequenceWithin(sequenceOut_)
3709            <<CoinMessageEol;
3710#ifdef COIN_DEVELOP
3711          printf("flag d\n");
3712#endif
3713          setFlagged(sequenceOut_);
3714          progress_.clearBadTimes();
3715         
3716          // Go to safe
3717          factorization_->pivotTolerance(0.99);
3718          forceFactorization_=1; // a bit drastic but ..
3719          type = 2;
3720          //assert (internalFactorize(1)==0);
3721          if (internalFactorize(1)) {
3722            CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
3723            CoinMemcpyN(savedSolution_+numberColumns_ ,
3724                        numberRows_,rowActivityWork_);
3725            CoinMemcpyN(savedSolution_ ,
3726                        numberColumns_,columnActivityWork_);
3727            // restore extra stuff
3728            int dummy;
3729            matrix_->generalExpanded(this,6,dummy);
3730            // debug
3731            int returnCode = internalFactorize(1);
3732            while (returnCode) {
3733              // ouch
3734              // switch off dense
3735              int saveDense = factorization_->denseThreshold();
3736              factorization_->setDenseThreshold(0);
3737              // Go to safe
3738              factorization_->pivotTolerance(0.99);
3739              // make sure will do safe factorization
3740              pivotVariable_[0]=-1;
3741              returnCode=internalFactorize(2);
3742              factorization_->setDenseThreshold(saveDense);
3743            }
3744            // get correct bounds on all variables
3745            resetFakeBounds();
3746          }
3747        }
3748      }
3749      if (problemStatus_!=-4||numberPivots>10)
3750        problemStatus_=-3;
3751    }
3752  } else {
3753    //printf("testing with accuracy of %g and status of %d\n",alphaAccuracy_,problemStatus_);
3754    //count_alpha++;
3755    //if ((count_alpha%5000)==0)
3756    //printf("count alpha %d\n",count_alpha);
3757  }
3758  // at this stage status is -3 or -4 if looks infeasible
3759  // get primal and dual solutions
3760#if 0
3761  {
3762    int numberTotal = numberRows_+numberColumns_;
3763    double * saveSol = CoinCopyOfArray(solution_,numberTotal);
3764    double * saveDj = CoinCopyOfArray(dj_,numberTotal);
3765    double tolerance = type ? 1.0e-4 : 1.0e-8;
3766    // always if values pass
3767    double saveObj=objectiveValue_;
3768    double sumPrimal = sumPrimalInfeasibilities_;
3769    int numberPrimal = numberPrimalInfeasibilities_;
3770    double sumDual = sumDualInfeasibilities_;
3771    int numberDual = numberDualInfeasibilities_;
3772    gutsOfSolution(givenDuals,NULL);
3773    int j;
3774    double largestPrimal = tolerance;
3775    int iPrimal=-1;
3776    for (j=0;j<numberTotal;j++) {
3777      double difference = solution_[j]-saveSol[j];
3778      if (fabs(difference)>largestPrimal) {
3779        iPrimal=j;
3780        largestPrimal=fabs(difference);
3781      }
3782    }
3783    double largestDual = tolerance;
3784    int iDual=-1;
3785    for (j=0;j<numberTotal;j++) {
3786      double difference = dj_[j]-saveDj[j];
3787      if (fabs(difference)>largestDual&&upper_[j]>lower_[j]) {
3788        iDual=j;
3789        largestDual=fabs(difference);
3790      }
3791    }
3792    if (!type) {
3793      if (fabs(saveObj-objectiveValue_)>1.0e-5||
3794          numberPrimal!=numberPrimalInfeasibilities_||numberPrimal!=1||
3795          fabs(sumPrimal-sumPrimalInfeasibilities_)>1.0e-5||iPrimal>=0||
3796          numberDual!=numberDualInfeasibilities_||numberDual!=0||
3797          fabs(sumDual-sumDualInfeasibilities_)>1.0e-5||iDual>=0)
3798        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",
3799               type,numberIterations_,numberPivots,
3800               numberPrimal,numberPrimalInfeasibilities_,sumPrimal,sumPrimalInfeasibilities_,
3801               largestPrimal,iPrimal,
3802               numberDual,numberDualInfeasibilities_,sumDual,sumDualInfeasibilities_,
3803               largestDual,iDual,
3804               saveObj,objectiveValue_);
3805    } else {
3806      if (fabs(saveObj-objectiveValue_)>1.0e-5||
3807          numberPrimalInfeasibilities_||iPrimal>=0||
3808          numberDualInfeasibilities_||iDual>=0)
3809        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",
3810               type,numberIterations_,numberPivots,
3811               numberPrimal,numberPrimalInfeasibilities_,sumPrimal,sumPrimalInfeasibilities_,
3812               largestPrimal,iPrimal,
3813               numberDual,numberDualInfeasibilities_,sumDual,sumDualInfeasibilities_,
3814               largestDual,iDual,
3815               saveObj,objectiveValue_);
3816    }
3817    delete [] saveSol;
3818    delete [] saveDj;
3819  }
3820#else
3821  if (type||ifValuesPass)
3822    gutsOfSolution(givenDuals,NULL);
3823#endif
3824  // If bad accuracy treat as singular
3825  if ((largestPrimalError_>1.0e15||largestDualError_>1.0e15)&&numberIterations_) {
3826    // restore previous basis
3827    unflagVariables = false;
3828    changeMade_++; // say something changed
3829    // Keep any flagged variables
3830    int i;
3831    for (i=0;i<numberRows_+numberColumns_;i++) {
3832      if (flagged(i))
3833        saveStatus_[i] |= 64; //say flagged
3834    }
3835    CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
3836    CoinMemcpyN(savedSolution_+numberColumns_ ,
3837           numberRows_,rowActivityWork_);
3838    CoinMemcpyN(savedSolution_ ,
3839           numberColumns_,columnActivityWork_);
3840    // restore extra stuff
3841    int dummy;
3842    matrix_->generalExpanded(this,6,dummy);
3843    // get correct bounds on all variables
3844    resetFakeBounds();
3845    // need to reject something
3846    char x = isColumn(sequenceOut_) ? 'C' :'R';
3847    handler_->message(CLP_SIMPLEX_FLAG,messages_)
3848      <<x<<sequenceWithin(sequenceOut_)
3849      <<CoinMessageEol;
3850#ifdef COIN_DEVELOP
3851    printf("flag e\n");
3852#endif
3853    setFlagged(sequenceOut_);
3854    progress_.clearBadTimes();
3855   
3856    // Go to safer
3857    double newTolerance = CoinMin(1.1*factorization_->pivotTolerance(),0.99);
3858    factorization_->pivotTolerance(newTolerance);
3859    forceFactorization_=1; // a bit drastic but ..
3860    if (alphaAccuracy_!=-1.0)
3861      alphaAccuracy_=-2.0;
3862    type = 2;
3863    //assert (internalFactorize(1)==0);
3864    if (internalFactorize(1)) {
3865      CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
3866      CoinMemcpyN(savedSolution_+numberColumns_ ,
3867              numberRows_,rowActivityWork_);
3868      CoinMemcpyN(savedSolution_ ,
3869           numberColumns_,columnActivityWork_);
3870      // restore extra stuff
3871      int dummy;
3872      matrix_->generalExpanded(this,6,dummy);
3873      // debug
3874      int returnCode = internalFactorize(1);
3875      while (returnCode) {
3876        // ouch
3877        // switch off dense
3878        int saveDense = factorization_->denseThreshold();
3879        factorization_->setDenseThreshold(0);
3880        // Go to safe
3881        factorization_->pivotTolerance(0.99);
3882        // make sure will do safe factorization
3883        pivotVariable_[0]=-1;
3884        returnCode=internalFactorize(2);
3885        factorization_->setDenseThreshold(saveDense);
3886      }
3887      // get correct bounds on all variables
3888      resetFakeBounds();
3889    }
3890    // get primal and dual solutions
3891    gutsOfSolution(givenDuals,NULL);
3892  } else if (goodAccuracy()) {
3893    // Can reduce tolerance
3894    double newTolerance = CoinMax(0.99*factorization_->pivotTolerance(),saveData.pivotTolerance_);
3895    factorization_->pivotTolerance(newTolerance);
3896  } 
3897  // Double check infeasibility if no action
3898  if (progress_.lastIterationNumber(0)==numberIterations_) {
3899    if (dualRowPivot_->looksOptimal()) {
3900      numberPrimalInfeasibilities_ = 0;
3901      sumPrimalInfeasibilities_ = 0.0;
3902    }
3903#if 1
3904  } else {
3905    double thisObj = objectiveValue_;
3906    double lastObj = progress_.lastObjective(0);
3907    if(!ifValuesPass&&firstFree_<0) {
3908      if (lastObj>thisObj+1.0e-3*CoinMax(fabs(thisObj),fabs(lastObj))+1.0) {
3909        int maxFactor = factorization_->maximumPivots();
3910        if (maxFactor>10&&numberPivots>1) {
3911          //printf("lastobj %g thisobj %g\n",lastObj,thisObj);
3912          //if (forceFactorization_<0)
3913          //forceFactorization_= maxFactor;
3914          //forceFactorization_ = CoinMax(1,(forceFactorization_>>1));
3915          forceFactorization_=1;
3916          //printf("Reducing factorization frequency - bad backwards\n");
3917          unflagVariables = false;
3918          changeMade_++; // say something changed
3919          CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
3920          CoinMemcpyN(savedSolution_+numberColumns_ ,
3921                      numberRows_,rowActivityWork_);
3922          CoinMemcpyN(savedSolution_ ,
3923                      numberColumns_,columnActivityWork_);
3924          // restore extra stuff
3925          int dummy;
3926          matrix_->generalExpanded(this,6,dummy);
3927          // get correct bounds on all variables
3928          resetFakeBounds();
3929          if(factorization_->pivotTolerance()<0.2)
3930            factorization_->pivotTolerance(0.2);
3931          if (alphaAccuracy_!=-1.0)
3932            alphaAccuracy_=-2.0;
3933          if (internalFactorize(1)) {
3934            CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
3935            CoinMemcpyN(savedSolution_+numberColumns_ ,
3936                        numberRows_,rowActivityWork_);
3937            CoinMemcpyN(savedSolution_ ,
3938                        numberColumns_,columnActivityWork_);
3939            // restore extra stuff
3940            int dummy;
3941            matrix_->generalExpanded(this,6,dummy);
3942            // debug
3943            int returnCode = internalFactorize(1);
3944            while (returnCode) {
3945              // ouch
3946              // switch off dense
3947              int saveDense = factorization_->denseThreshold();
3948              factorization_->setDenseThreshold(0);
3949              // Go to safe
3950              factorization_->pivotTolerance(0.99);
3951              // make sure will do safe factorization
3952              pivotVariable_[0]=-1;
3953              returnCode=internalFactorize(2);
3954              factorization_->setDenseThreshold(saveDense);
3955            }
3956            resetFakeBounds();
3957          }
3958          type = 2; // so will restore weights
3959          // get primal and dual solutions
3960          gutsOfSolution(givenDuals,NULL);
3961        }
3962      } else if (lastObj<thisObj-1.0e-5*CoinMax(fabs(thisObj),fabs(lastObj))-1.0e-3) {
3963        numberTimesOptimal_=0;
3964      }
3965    }
3966#endif
3967  }
3968  // Up tolerance if looks a bit odd
3969  if (numberIterations_>CoinMax(1000,numberRows_>>4)&&(specialOptions_&64)!=0) {
3970    if (sumPrimalInfeasibilities_&&sumPrimalInfeasibilities_<1.0e5) {
3971      int backIteration = progress_.lastIterationNumber(CLP_PROGRESS-1);
3972      if (backIteration>0&&numberIterations_-backIteration<9*CLP_PROGRESS) {
3973        if (factorization_->pivotTolerance()<0.9) {
3974          // up tolerance
3975          factorization_->pivotTolerance(CoinMin(factorization_->pivotTolerance()*1.05+0.02,0.91));
3976          //printf("tol now %g\n",factorization_->pivotTolerance());
3977          progress_.clearIterationNumbers();
3978        }
3979      }
3980    }
3981  }
3982  // Check if looping
3983  int loop;
3984  if (!givenDuals&&type!=2) 
3985    loop = progress_.looping();
3986  else
3987    loop=-1;
3988  int situationChanged=0;
3989  if (loop>=0) {
3990    problemStatus_ = loop; //exit if in loop
3991    if (!problemStatus_) {
3992      // declaring victory
3993      numberPrimalInfeasibilities_ = 0;
3994      sumPrimalInfeasibilities_ = 0.0;
3995    } else {
3996      problemStatus_ = 10; // instead - try other algorithm
3997#if COIN_DEVELOP>2
3998      printf("returning at %d\n",__LINE__);
3999#endif
4000    }
4001    return;
4002  } else if (loop<-1) {
4003    // something may have changed
4004    gutsOfSolution(NULL,NULL);
4005    situationChanged=1;
4006  }
4007  // really for free variables in
4008  if((progressFlag_&2)!=0) {
4009    situationChanged=2;
4010  }
4011  progressFlag_ = 0; //reset progress flag
4012  /*if (!numberIterations_&&sumDualInfeasibilities_)
4013    printf("OBJ %g sumPinf %g sumDinf %g\n",
4014           objectiveValue(),sumPrimalInfeasibilities_,
4015           sumDualInfeasibilities_);*/
4016  if (handler_->detail(CLP_SIMPLEX_STATUS,messages_)<100) {
4017    handler_->message(CLP_SIMPLEX_STATUS,messages_)
4018      <<numberIterations_<<objectiveValue();
4019    handler_->printing(sumPrimalInfeasibilities_>0.0)
4020      <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
4021    handler_->printing(sumDualInfeasibilities_>0.0)
4022      <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
4023    handler_->printing(numberDualInfeasibilitiesWithoutFree_
4024                       <numberDualInfeasibilities_)
4025                         <<numberDualInfeasibilitiesWithoutFree_;
4026    handler_->message()<<CoinMessageEol;
4027  }
4028  realDualInfeasibilities=sumDualInfeasibilities_;
4029  double saveTolerance =dualTolerance_;
4030  // If we need to carry on cleaning variables
4031  if (!numberPrimalInfeasibilities_&&(specialOptions_&1024)!=0&&CLEAN_FIXED) {
4032    for (int iRow=0;iRow<numberRows_;iRow++) {
4033      int iPivot=pivotVariable_[iRow];
4034      if (!flagged(iPivot)&&pivoted(iPivot)) {
4035        // carry on
4036        numberPrimalInfeasibilities_=-1;
4037        sumOfRelaxedPrimalInfeasibilities_ = 1.0;
4038        sumPrimalInfeasibilities_ = 1.0;
4039        break;
4040      }
4041    }
4042  }
4043  /* If we are primal feasible and any dual infeasibilities are on
4044     free variables then it is better to go to primal */
4045  if (!numberPrimalInfeasibilities_&&!numberDualInfeasibilitiesWithoutFree_&&
4046      numberDualInfeasibilities_)
4047    problemStatus_=10;
4048  // dual bound coming in
4049  //double saveDualBound = dualBound_;
4050  while (problemStatus_<=-3) {
4051    int cleanDuals=0;
4052    if (situationChanged!=0)
4053      cleanDuals=1;
4054    int numberChangedBounds=0;
4055    int doOriginalTolerance=0;
4056    if ( lastCleaned==numberIterations_)
4057      doOriginalTolerance=1;
4058    // check optimal
4059    // give code benefit of doubt
4060    if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
4061        sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
4062      // say optimal (with these bounds etc)
4063      numberDualInfeasibilities_ = 0;
4064      sumDualInfeasibilities_ = 0.0;
4065      numberPrimalInfeasibilities_ = 0;
4066      sumPrimalInfeasibilities_ = 0.0;
4067    }
4068    //if (dualFeasible()||problemStatus_==-4||(primalFeasible()&&!numberDualInfeasibilitiesWithoutFree_)) {
4069    if (dualFeasible()||problemStatus_==-4) {
4070      progress_.modifyObjective(objectiveValue_
4071                               -sumDualInfeasibilities_*dualBound_);
4072      // see if cutoff reached
4073      double limit = 0.0;
4074      getDblParam(ClpDualObjectiveLimit, limit);
4075      if(fabs(limit)<1.0e30&&objectiveValue()*optimizationDirection_>
4076         limit+1.0e-7+1.0e-8*fabs(limit)&&!numberAtFakeBound()) {
4077        //looks infeasible on objective
4078        if (perturbation_==101) {
4079          perturbation_=102; // stop any perturbations
4080          cleanDuals=1;
4081          // make sure fake bounds are back
4082          changeBounds(1,NULL,changeCost);
4083          createRim4(false);
4084          // make sure duals are current
4085          computeDuals(givenDuals);
4086          checkDualSolution();
4087        }
4088      }
4089      if (primalFeasible()&&!givenDuals) {
4090        // may be optimal - or may be bounds are wrong
4091        handler_->message(CLP_DUAL_BOUNDS,messages_)
4092          <<dualBound_
4093          <<CoinMessageEol;
4094        // save solution in case unbounded
4095        double * saveColumnSolution = NULL;
4096        double * saveRowSolution = NULL;
4097        bool inCbc = (specialOptions_&(0x01000000|16384))!=0;
4098        if (!inCbc) {
4099          saveColumnSolution = CoinCopyOfArray(columnActivityWork_,numberColumns_);
4100          saveRowSolution = CoinCopyOfArray(rowActivityWork_,numberRows_);
4101        }
4102        numberChangedBounds=changeBounds(0,rowArray_[3],changeCost);
4103        if (numberChangedBounds<=0&&!numberDualInfeasibilities_) {
4104          //looks optimal - do we need to reset tolerance
4105          if (perturbation_==101) {
4106            perturbation_=102; // stop any perturbations
4107            cleanDuals=1;
4108            // make sure fake bounds are back
4109            //computeObjectiveValue();
4110            changeBounds(1,NULL,changeCost);
4111            //computeObjectiveValue();
4112            createRim4(false);
4113            // make sure duals are current
4114            computeDuals(givenDuals);
4115            checkDualSolution();
4116#define DUAL_TRY_FASTER
4117#ifdef DUAL_TRY_FASTER
4118            if (numberDualInfeasibilities_) {
4119#endif
4120              numberChanged_=1; // force something to happen
4121              lastCleaned=numberIterations_-1;
4122#ifdef DUAL_TRY_FASTER
4123            } else {
4124              //double value = objectiveValue_;
4125              computeObjectiveValue(true);
4126              //printf("old %g new %g\n",value,objectiveValue_);
4127              //numberChanged_=1;
4128            }
4129#endif
4130          }
4131          if (lastCleaned<numberIterations_&&numberTimesOptimal_<4&&
4132              (numberChanged_||(specialOptions_&4096)==0)) {
4133            doOriginalTolerance=2;
4134            numberTimesOptimal_++;
4135            changeMade_++; // say something changed
4136            if (numberTimesOptimal_==1) {
4137              dualTolerance_ = dblParam_[ClpDualTolerance];
4138            } else {
4139              if (numberTimesOptimal_==2) {
4140                // better to have small tolerance even if slower
4141                factorization_->zeroTolerance(1.0e-15);
4142              }
4143              dualTolerance_ = dblParam_[ClpDualTolerance];
4144              dualTolerance_ *= pow(2.0,numberTimesOptimal_-1);
4145            }
4146            cleanDuals=2; // If nothing changed optimal else primal
4147          } else {
4148            problemStatus_=0; // optimal
4149            if (lastCleaned<numberIterations_&&numberChanged_) {
4150              handler_->message(CLP_SIMPLEX_GIVINGUP,messages_)
4151                <<CoinMessageEol;
4152            }
4153          }
4154        } else {
4155          cleanDuals=1;
4156          if (doOriginalTolerance==1) {
4157            // check unbounded
4158            // find a variable with bad dj
4159            int iSequence;
4160            int iChosen=-1;
4161            if (!inCbc) {
4162              double largest = 100.0*primalTolerance_;
4163              for (iSequence=0;iSequence<numberRows_+numberColumns_;
4164                   iSequence++) {
4165                double djValue = dj_[iSequence];
4166                double originalLo = originalLower(iSequence);
4167                double originalUp = originalUpper(iSequence);
4168                if (fabs(djValue)>fabs(largest)) {
4169                  if (getStatus(iSequence)!=basic) {
4170                    if (djValue>0&&originalLo<-1.0e20) {
4171                      if (djValue>fabs(largest)) {
4172                        largest=djValue;
4173                        iChosen=iSequence;
4174                      }
4175                    } else if (djValue<0&&originalUp>1.0e20) {
4176                      if (-djValue>fabs(largest)) {
4177                        largest=djValue;
4178                        iChosen=iSequence;
4179                      }
4180                    } 
4181                  }
4182                }
4183              }
4184            }
4185            if (iChosen>=0) {
4186              int iSave=sequenceIn_;
4187              sequenceIn_=iChosen;
4188              unpack(rowArray_[1]);
4189              sequenceIn_ = iSave;
4190              // if dual infeasibilities then must be free vector so add in dual
4191              if (numberDualInfeasibilities_) {
4192                if (fabs(changeCost)>1.0e-5)
4193                  printf("Odd free/unbounded combo\n");
4194                changeCost += cost_[iChosen];
4195              }
4196              problemStatus_ = checkUnbounded(rowArray_[1],rowArray_[0],
4197                                              changeCost);
4198              rowArray_[1]->clear();
4199            } else {
4200              problemStatus_=-3;
4201            }
4202            if (problemStatus_==2&&perturbation_==101) {
4203              perturbation_=102; // stop any perturbations
4204              cleanDuals=1;
4205              createRim4(false);
4206              problemStatus_=-1;
4207            }
4208            if (problemStatus_==2) {
4209              // it is unbounded - restore solution
4210              // but first add in changes to non-basic
4211              int iColumn;
4212              double * original = columnArray_[0]->denseVector();
4213              for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4214                if(getColumnStatus(iColumn)!= basic)
4215                  ray_[iColumn] += 
4216                    saveColumnSolution[iColumn]-original[iColumn];
4217                columnActivityWork_[iColumn] = original[iColumn];
4218              }
4219              CoinMemcpyN(saveRowSolution,numberRows_,
4220                                rowActivityWork_);
4221            }
4222          } else {
4223            doOriginalTolerance=2;
4224            rowArray_[0]->clear();
4225          }
4226        }
4227        delete [] saveColumnSolution;
4228        delete [] saveRowSolution;
4229      }
4230      if (problemStatus_==-4||problemStatus_==-5) {
4231        // may be infeasible - or may be bounds are wrong
4232        numberChangedBounds=changeBounds(0,NULL,changeCost);
4233        /* Should this be here as makes no difference to being feasible.
4234           But seems to make a difference to run times. */
4235        if (perturbation_==101&&0) {
4236          perturbation_=102; // stop any perturbations
4237          cleanDuals=1;
4238          numberChangedBounds=1;
4239          // make sure fake bounds are back
4240          changeBounds(1,NULL,changeCost);
4241          createRim4(false);
4242        }
4243        if ((numberChangedBounds<=0||dualBound_>1.0e20||
4244            (largestPrimalError_>1.0&&dualBound_>1.0e17))&&
4245            (numberPivots<4||sumPrimalInfeasibilities_>1.0e-6)) {
4246          problemStatus_=1; // infeasible
4247          if (perturbation_==101) {
4248            perturbation_=102; // stop any perturbations
4249            //cleanDuals=1;
4250            //numberChangedBounds=1;
4251            //createRim4(false);
4252          }
4253        } else {
4254          problemStatus_=-1; //iterate
4255          cleanDuals=1;
4256          if (numberChangedBounds<=0)
4257            doOriginalTolerance=2;
4258          // and delete ray which has been created
4259          delete [] ray_;
4260          ray_ = NULL;
4261        }
4262
4263      }
4264    } else {
4265      cleanDuals=1;
4266    }
4267    if (problemStatus_<0) {
4268      if (doOriginalTolerance==2) {
4269        // put back original tolerance
4270        lastCleaned=numberIterations_;
4271        numberChanged_ =0; // Number of variables with changed costs
4272        handler_->message(CLP_DUAL_ORIGINAL,messages_)
4273          <<CoinMessageEol;
4274        perturbation_=102; // stop any perturbations
4275#if 0
4276        double * xcost = new double[numberRows_+numberColumns_];
4277        double * xlower = new double[numberRows_+numberColumns_];
4278        double * xupper = new double[numberRows_+numberColumns_];
4279        double * xdj = new double[numberRows_+numberColumns_];
4280        double * xsolution = new double[numberRows_+numberColumns_];
4281 CoinMemcpyN(cost_,(numberRows_+numberColumns_),xcost);
4282 CoinMemcpyN(lower_,(numberRows_+numberColumns_),xlower);
4283 CoinMemcpyN(upper_,(numberRows_+numberColumns_),xupper);
4284 CoinMemcpyN(dj_,(numberRows_+numberColumns_),xdj);
4285 CoinMemcpyN(solution_,(numberRows_+numberColumns_),xsolution);
4286#endif
4287        createRim4(false);
4288        // make sure duals are current
4289        computeDuals(givenDuals);
4290        checkDualSolution();
4291#if 0
4292        int i;
4293        for (i=0;i<numberRows_+numberColumns_;i++) {
4294          if (cost_[i]!=xcost[i])
4295            printf("** %d old cost %g new %g sol %g\n",
4296                   i,xcost[i],cost_[i],solution_[i]);
4297          if (lower_[i]!=xlower[i])
4298            printf("** %d old lower %g new %g sol %g\n",
4299                   i,xlower[i],lower_[i],solution_[i]);
4300          if (upper_[i]!=xupper[i])
4301            printf("** %d old upper %g new %g sol %g\n",
4302                   i,xupper[i],upper_[i],solution_[i]);
4303          if (dj_[i]!=xdj[i])
4304            printf("** %d old dj %g new %g sol %g\n",
4305                   i,xdj[i],dj_[i],solution_[i]);
4306          if (solution_[i]!=xsolution[i])
4307            printf("** %d old solution %g new %g sol %g\n",
4308                   i,xsolution[i],solution_[i],solution_[i]);
4309        }
4310        //delete [] xcost;
4311        //delete [] xupper;
4312        //delete [] xlower;
4313        //delete [] xdj;
4314        //delete [] xsolution;
4315#endif
4316        // put back bounds as they were if was optimal
4317        if (doOriginalTolerance==2&&cleanDuals!=2) {
4318          changeMade_++; // say something changed
4319          /* We may have already changed some bounds in this function
4320             so save numberFake_ and add in.
4321
4322             Worst that can happen is that we waste a bit of time  - but it must be finite.
4323          */
4324          int saveNumberFake = numberFake_;
4325          changeBounds(1,NULL,changeCost);
4326          numberFake_ += saveNumberFake;
4327          cleanDuals=2;
4328          //cleanDuals=1;
4329        }
4330#if 0
4331        //int i;
4332        for (i=0;i<numberRows_+numberColumns_;i++) {
4333          if (cost_[i]!=xcost[i])
4334            printf("** %d old cost %g new %g sol %g\n",
4335                   i,xcost[i],cost_[i],solution_[i]);
4336          if (lower_[i]!=xlower[i])
4337            printf("** %d old lower %g new %g sol %g\n",
4338                   i,xlower[i],lower_[i],solution_[i]);
4339          if (upper_[i]!=xupper[i])
4340            printf("** %d old upper %g new %g sol %g\n",
4341                   i,xupper[i],upper_[i],solution_[i]);
4342          if (dj_[i]!=xdj[i])
4343            printf("** %d old dj %g new %g sol %g\n",
4344                   i,xdj[i],dj_[i],solution_[i]);
4345          if (solution_[i]!=xsolution[i])
4346            printf("** %d old solution %g new %g sol %g\n",
4347                   i,xsolution[i],solution_[i],solution_[i]);
4348        }
4349        delete [] xcost;
4350        delete [] xupper;
4351        delete [] xlower;
4352        delete [] xdj;
4353        delete [] xsolution;
4354#endif
4355      }
4356      if (cleanDuals==1||(cleanDuals==2&&!numberDualInfeasibilities_)) {
4357        // make sure dual feasible
4358        // look at all rows and columns
4359        rowArray_[0]->clear();
4360        columnArray_[0]->clear();
4361        double objectiveChange=0.0;
4362#if 0
4363        double * xcost = new double[numberRows_+numberColumns_];
4364        double * xlower = new double[numberRows_+numberColumns_];
4365        double * xupper = new double[numberRows_+numberColumns_];
4366        double * xdj = new double[numberRows_+numberColumns_];
4367        double * xsolution = new double[numberRows_+numberColumns_];
4368 CoinMemcpyN(cost_,(numberRows_+numberColumns_),xcost);
4369 CoinMemcpyN(lower_,(numberRows_+numberColumns_),xlower);
4370 CoinMemcpyN(upper_,(numberRows_+numberColumns_),xupper);
4371 CoinMemcpyN(dj_,(numberRows_+numberColumns_),xdj);
4372 CoinMemcpyN(solution_,(numberRows_+numberColumns_),xsolution);
4373#endif
4374        if (givenDuals)
4375          dualTolerance_=1.0e50;
4376        updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[1],
4377          0.0,objectiveChange,true);
4378        dualTolerance_=saveTolerance;
4379#if 0
4380        int i;
4381        for (i=0;i<numberRows_+numberColumns_;i++) {
4382          if (cost_[i]!=xcost[i])
4383            printf("** %d old cost %g new %g sol %g\n",
4384                   i,xcost[i],cost_[i],solution_[i]);
4385          if (lower_[i]!=xlower[i])
4386            printf("** %d old lower %g new %g sol %g\n",
4387                   i,xlower[i],lower_[i],solution_[i]);
4388          if (upper_[i]!=xupper[i])
4389            printf("** %d old upper %g new %g sol %g\n",
4390                   i,xupper[i],upper_[i],solution_[i]);
4391          if (dj_[i]!=xdj[i])
4392            printf("** %d old dj %g new %g sol %g\n",
4393                   i,xdj[i],dj_[i],solution_[i]);
4394          if (solution_[i]!=xsolution[i])
4395            printf("** %d old solution %g new %g sol %g\n",
4396                   i,xsolution[i],solution_[i],solution_[i]);
4397        }
4398        delete [] xcost;
4399        delete [] xupper;
4400        delete [] xlower;
4401        delete [] xdj;
4402        delete [] xsolution;
4403#endif
4404        // for now - recompute all
4405        gutsOfSolution(NULL,NULL);
4406        if (givenDuals)
4407          dualTolerance_=1.0e50;
4408        updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[1],
4409          0.0,objectiveChange,true);
4410        dualTolerance_=saveTolerance;
4411        //assert(numberDualInfeasibilitiesWithoutFree_==0);
4412        if (numberDualInfeasibilities_) {
4413          if (numberPrimalInfeasibilities_||numberPivots) {
4414            problemStatus_=-1; // carry on as normal
4415          } else {
4416            problemStatus_=10; // try primal
4417#if COIN_DEVELOP>1
4418            printf("returning at %d\n",__LINE__);
4419#endif
4420          }
4421        } else if (situationChanged==2) {
4422          problemStatus_=-1; // carry on as normal
4423        }
4424        situationChanged=0;
4425      } else {
4426        // iterate
4427        if (cleanDuals!=2) {
4428          problemStatus_=-1;
4429        } else {
4430          problemStatus_ = 10; // try primal
4431#if COIN_DEVELOP>2
4432          printf("returning at %d\n",__LINE__);
4433#endif
4434        }
4435      }
4436    }
4437  }
4438  // unflag all variables (we may want to wait a bit?)
4439  if ((tentativeStatus!=-2&&tentativeStatus!=-1)&&unflagVariables) {
4440    int iRow;
4441    int numberFlagged=0;
4442    for (iRow=0;iRow<numberRows_;iRow++) {
4443      int iPivot=pivotVariable_[iRow];
4444      if (flagged(iPivot)) {
4445        numberFlagged++;
4446        clearFlagged(iPivot);
4447      }
4448    }
4449#ifdef COIN_DEVELOP
4450    if (numberFlagged) {
4451      printf("unflagging %d variables - tentativeStatus %d probStat %d ninf %d nopt %d\n",numberFlagged,tentativeStatus,
4452             problemStatus_,numberPrimalInfeasibilities_,
4453             numberTimesOptimal_);
4454    }
4455#endif
4456    unflagVariables = numberFlagged>0;
4457    if (numberFlagged&&!numberPivots) {
4458      /* looks like trouble as we have not done any iterations.
4459         Try changing pivot tolerance then give it a few goes and give up */
4460      if (factorization_->pivotTolerance()<0.9) {
4461        factorization_->pivotTolerance(0.99);
4462        problemStatus_=-1;
4463      } else if (numberTimesOptimal_<3) {
4464        numberTimesOptimal_++;
4465        problemStatus_=-1;
4466      } else {
4467        unflagVariables=false;
4468        //secondaryStatus_ = 1; // and say probably infeasible
4469        // try primal
4470        problemStatus_=10;
4471#if COIN_DEVELOP>1
4472        printf("returning at %d\n",__LINE__);
4473#endif
4474      }
4475    }
4476  }
4477  if (problemStatus_<0) {
4478    if (type==0||type==1) {
4479      if (!type) {
4480        // create save arrays
4481        delete [] saveStatus_;
4482        delete [] savedSolution_;
4483        saveStatus_ = new unsigned char [numberRows_+numberColumns_];
4484        savedSolution_ = new double [numberRows_+numberColumns_];
4485      }
4486      // save arrays
4487      CoinMemcpyN(status_,numberColumns_+numberRows_,saveStatus_);
4488      CoinMemcpyN(rowActivityWork_,
4489                  numberRows_,savedSolution_+numberColumns_);
4490      CoinMemcpyN(columnActivityWork_,numberColumns_,savedSolution_);
4491      // save extra stuff
4492      int dummy;
4493      matrix_->generalExpanded(this,5,dummy);
4494    }
4495    if (weightsSaved) {
4496      // restore weights (if saved) - also recompute infeasibility list
4497      if (tentativeStatus>-3) 
4498        dualRowPivot_->saveWeights(this,(type <2) ? 2 : 4);
4499      else
4500        dualRowPivot_->saveWeights(this,3);
4501    }
4502  }
4503  // see if cutoff reached
4504  double limit = 0.0;
4505  getDblParam(ClpDualObjectiveLimit, limit);
4506#if 0
4507  if(fabs(limit)<1.0e30&&objectiveValue()*optimizationDirection_>
4508     limit+100.0) {
4509    printf("lim %g obj %g %g - wo perturb %g sum dual %g\n",
4510           limit,objectiveValue_,objectiveValue(),computeInternalObjectiveValue(),sumDualInfeasibilities_);
4511  }
4512#endif
4513  if(fabs(limit)<1.0e30&&objectiveValue()*optimizationDirection_>
4514     limit&&!numberAtFakeBound()) {
4515    bool looksInfeasible = !numberDualInfeasibilities_;
4516    if (objectiveValue()*optimizationDirection_>limit+fabs(0.1*limit)+1.0e2*sumDualInfeasibilities_+1.0e4&&
4517        sumDualInfeasibilities_<largestDualError_&&numberIterations_>0.5*numberRows_+1000)
4518      looksInfeasible=true;
4519    if (looksInfeasible) {
4520      // Even if not perturbed internal costs may have changed
4521      // be careful
4522      if (numberIterations_) {
4523        if(computeInternalObjectiveValue()>limit) {
4524          problemStatus_=1;
4525          secondaryStatus_ = 1; // and say was on cutoff
4526        }
4527      } else {
4528        problemStatus_=1;
4529        secondaryStatus_ = 1; // and say was on cutoff
4530      }
4531    }
4532  }
4533  // If we are in trouble and in branch and bound give up
4534  if ((specialOptions_&1024)!=0) {
4535    int looksBad=0;
4536    if (largestPrimalError_*largestDualError_>1.0e2) {
4537      looksBad=1;
4538    } else if (largestPrimalError_>1.0e-2
4539        &&objectiveValue_>CoinMin(1.0e15,1.0e3*limit)) {
4540      looksBad=2;
4541    }
4542    if (looksBad) {
4543      if (factorization_->pivotTolerance()<0.9) {
4544        // up tolerance
4545        factorization_->pivotTolerance(CoinMin(factorization_->pivotTolerance()*1.05+0.02,0.91));
4546      } else if (numberIterations_>10000) {
4547        if (handler_->logLevel()>2)
4548          printf("bad dual - saying infeasible %d\n",looksBad);
4549        problemStatus_=1;
4550        secondaryStatus_ = 1; // and say was on cutoff
4551      } else if (largestPrimalError_>1.0e5) {
4552        {
4553          int iBigB=-1;
4554          double bigB=0.0;
4555          int iBigN=-1;
4556          double bigN=0.0;
4557          for (int i=0;i<numberRows_+numberColumns_;i++) {
4558            double value = fabs(solution_[i]);
4559            if (getStatus(i)==basic) {
4560              if (value>bigB) {
4561                bigB= value;
4562                iBigB=i;
4563              }
4564            } else {
4565              if (value>bigN) {
4566                bigN= value;
4567                iBigN=i;
4568              }
4569            }
4570          }
4571#ifdef CLP_INVESTIGATE
4572          if (bigB>1.0e8||bigN>1.0e8) {
4573            if (handler_->logLevel()>0)
4574              printf("it %d - basic %d %g, nonbasic %d %g\n",
4575                     numberIterations_,iBigB,bigB,iBigN,bigN);
4576          }
4577#endif
4578        }
4579#if COIN_DEVELOP!=2
4580        if (handler_->logLevel()>2)
4581#endif
4582          printf("bad dual - going to primal %d %g\n",looksBad,largestPrimalError_);
4583        allSlackBasis(true);
4584        problemStatus_=10;
4585      }
4586    }
4587  }
4588  if (problemStatus_<0&&!changeMade_) {
4589    problemStatus_=4; // unknown
4590  }
4591  lastGoodIteration_ = numberIterations_;
4592  if (numberIterations_>lastBadIteration_+100)
4593    moreSpecialOptions_ &= ~16; // clear check accuracy flag
4594  if (problemStatus_<0) {
4595    sumDualInfeasibilities_=realDualInfeasibilities; // back to say be careful
4596    if (sumDualInfeasibilities_)
4597      numberDualInfeasibilities_=1;
4598  }
4599#if 1
4600  double thisObj = progress_.lastObjective(0);
4601  double lastObj = progress_.lastObjective(1);
4602  if (lastObj>thisObj+1.0e-4*CoinMax(fabs(thisObj),fabs(lastObj))+1.0e-4
4603      &&givenDuals==NULL&&firstFree_<0) {
4604    int maxFactor = factorization_->maximumPivots();
4605    if (maxFactor>10) {
4606      if (forceFactorization_<0)
4607        forceFactorization_= maxFactor;
4608      forceFactorization_ = CoinMax(1,(forceFactorization_>>1));
4609      //printf("Reducing factorization frequency\n");
4610    } 
4611  }
4612#endif
4613  // Allow matrices to be sorted etc
4614  int fake=-999; // signal sort
4615  matrix_->correctSequence(this,fake,fake);
4616  if (alphaAccuracy_>0.0)
4617      alphaAccuracy_=1.0;
4618}
4619/* While updateDualsInDual sees what effect is of flip
4620   this does actual flipping.
4621   If change >0.0 then value in array >0.0 => from lower to upper
4622*/
4623void 
4624ClpSimplexDual::flipBounds(CoinIndexedVector * rowArray,
4625                  CoinIndexedVector * columnArray,
4626                  double change)
4627{
4628  int number;
4629  int * which;
4630 
4631  int iSection;
4632
4633  for (iSection=0;iSection<2;iSection++) {
4634    int i;
4635    double * solution = solutionRegion(iSection);
4636    double * lower = lowerRegion(iSection);
4637    double * upper = upperRegion(iSection);
4638    int addSequence;
4639    if (!iSection) {
4640      number = rowArray->getNumElements();
4641      which = rowArray->getIndices();
4642      addSequence = numberColumns_;
4643    } else {
4644      number = columnArray->getNumElements();
4645      which = columnArray->getIndices();
4646      addSequence = 0;
4647    }
4648   
4649    for (i=0;i<number;i++) {
4650      int iSequence = which[i];
4651      Status status = getStatus(iSequence+addSequence);
4652
4653      switch(status) {
4654
4655      case basic:
4656      case isFree:
4657      case superBasic:
4658      case ClpSimplex::isFixed:
4659        break;
4660      case atUpperBound:
4661        // to lower bound
4662        setStatus(iSequence+addSequence,atLowerBound);
4663        solution[iSequence] = lower[iSequence];
4664        break;
4665      case atLowerBound:
4666        // to upper bound
4667        setStatus(iSequence+addSequence,atUpperBound);
4668        solution[iSequence] = upper[iSequence];
4669        break;
4670      }
4671    }
4672  }
4673  rowArray->setNumElements(0);
4674  columnArray->setNumElements(0);
4675}
4676// Restores bound to original bound
4677void 
4678ClpSimplexDual::originalBound( int iSequence)
4679{
4680  if (getFakeBound(iSequence)!=noFake) {
4681    numberFake_--;;
4682    setFakeBound(iSequence,noFake);
4683#ifdef CLP_AUXILIARY_MODEL
4684    if (auxiliaryModel_) {
4685      // just copy back
4686      lower_[iSequence]=auxiliaryModel_->lowerRegion()[iSequence+numberRows_+numberColumns_];
4687      upper_[iSequence]=auxiliaryModel_->upperRegion()[iSequence+numberRows_+numberColumns_];
4688      return;
4689    }
4690#endif
4691    if (iSequence>=numberColumns_) {
4692      // rows
4693      int iRow = iSequence-numberColumns_;
4694      rowLowerWork_[iRow]=rowLower_[iRow];
4695      rowUpperWork_[iRow]=rowUpper_[iRow];
4696      if (rowScale_) {
4697        if (rowLowerWork_[iRow]>-1.0e50)
4698          rowLowerWork_[iRow] *= rowScale_[iRow]*rhsScale_;
4699        if (rowUpperWork_[iRow]<1.0e50)
4700          rowUpperWork_[iRow] *= rowScale_[iRow]*rhsScale_;
4701      } else if (rhsScale_!=1.0) {
4702        if (rowLowerWork_[iRow]>-1.0e50)
4703          rowLowerWork_[iRow] *= rhsScale_;
4704        if (rowUpperWork_[iRow]<1.0e50)
4705          rowUpperWork_[iRow] *= rhsScale_;
4706      }
4707    } else {
4708      // columns
4709      columnLowerWork_[iSequence]=columnLower_[iSequence];
4710      columnUpperWork_[iSequence]=columnUpper_[iSequence];
4711      if (rowScale_) {
4712        double multiplier = 1.0/columnScale_[iSequence];
4713        if (columnLowerWork_[iSequence]>-1.0e50)
4714          columnLowerWork_[iSequence] *= multiplier*rhsScale_;
4715        if (columnUpperWork_[iSequence]<1.0e50)
4716          columnUpperWork_[iSequence] *= multiplier*rhsScale_;
4717      } else if (rhsScale_!=1.0) {
4718        if (columnLowerWork_[iSequence]>-1.0e50)
4719          columnLowerWork_[iSequence] *= rhsScale_;
4720        if (columnUpperWork_[iSequence]<1.0e50)
4721          columnUpperWork_[iSequence] *= rhsScale_;
4722      }
4723    }
4724  }
4725}
4726/* As changeBounds but just changes new bounds for a single variable.
4727   Returns true if change */
4728bool 
4729ClpSimplexDual::changeBound( int iSequence)
4730{
4731  // old values
4732  double oldLower=lower_[iSequence];
4733  double oldUpper=upper_[iSequence];
4734  double value=solution_[iSequence];
4735  bool modified=false;
4736  originalBound(iSequence);
4737  // original values
4738  double lowerValue=lower_[iSequence];
4739  double upperValue=upper_[iSequence];
4740  // back to altered values
4741  lower_[iSequence] = oldLower;
4742  upper_[iSequence] = oldUpper;
4743  assert (getFakeBound(iSequence)==noFake);
4744  //if (getFakeBound(iSequence)!=noFake)
4745  //numberFake_--;;
4746  if (value==oldLower) {
4747    if (upperValue > oldLower + dualBound_) {
4748      upper_[iSequence]=oldLower+dualBound_;
4749      setFakeBound(iSequence,upperFake);
4750      modified=true;
4751      numberFake_++;
4752    }
4753  } else if (value==oldUpper) {
4754    if (lowerValue < oldUpper - dualBound_) {
4755      lower_[iSequence]=oldUpper-dualBound_;
4756      setFakeBound(iSequence,lowerFake);
4757      modified=true;
4758      numberFake_++;
4759    }
4760  } else {
4761    assert(value==oldLower||value==oldUpper);
4762  }
4763  return modified;
4764}
4765// Perturbs problem
4766int 
4767ClpSimplexDual::perturb()
4768{
4769  if (perturbation_>100)
4770    return 0; //perturbed already
4771  if (perturbation_==100)
4772    perturbation_=50; // treat as normal
4773  int savePerturbation = perturbation_;
4774  bool modifyRowCosts=false;
4775  // dual perturbation
4776  double perturbation=1.0e-20;
4777  // maximum fraction of cost to perturb
4778  double maximumFraction = 1.0e-5;
4779  double constantPerturbation = 100.0*dualTolerance_;
4780  int maxLength=0;
4781  int minLength=numberRows_;
4782  double averageCost = 0.0;
4783#if 0
4784  // look at element range
4785  double smallestNegative;
4786  double largestNegative;
4787  double smallestPositive;
4788  double largestPositive;
4789  matrix_->rangeOfElements(smallestNegative, largestNegative,
4790                           smallestPositive, largestPositive);
4791  smallestPositive = CoinMin(fabs(smallestNegative),smallestPositive);
4792  largestPositive = CoinMax(fabs(largestNegative),largestPositive);
4793  double elementRatio = largestPositive/smallestPositive;
4794#endif
4795  int numberNonZero=0;
4796  if (!numberIterations_&&perturbation_>=50) {
4797    // See if we need to perturb
4798    double * sort = new double[numberColumns_];
4799    // Use objective BEFORE scaling
4800    const double * obj = objective();
4801    int i;
4802    for (i=0;i<numberColumns_;i++) {
4803      double value = fabs(obj[i]);
4804      sort[i]=value;
4805      averageCost += value;
4806      if (value)
4807        numberNonZero++;
4808    }
4809    if (numberNonZero)
4810      averageCost /= (double) numberNonZero;
4811    else
4812      averageCost = 1.0;
4813    std::sort(sort,sort+numberColumns_);
4814    int number=1;
4815    double last = sort[0];
4816    for (i=1;i<numberColumns_;i++) {
4817      if (last!=sort[i])
4818        number++;
4819      last=sort[i];
4820    }
4821    delete [] sort;
4822    if (!numberNonZero)
4823      return 1; // safer to use primal
4824#if 0
4825    printf("nnz %d percent %d",number,(number*100)/numberColumns_);
4826    if (number*4>numberColumns_)
4827      printf(" - Would not perturb\n");
4828    else
4829      printf(" - Would perturb\n");
4830    //exit(0);
4831#endif
4832    //printf("ratio number diff costs %g, element ratio %g\n",((double)number)/((double) numberColumns_),
4833    //                                                                elementRatio);
4834    //number=0;
4835    //if (number*4>numberColumns_||elementRatio>1.0e12) {
4836    if (number*4>numberColumns_) {
4837      perturbation_=100;
4838      return 0; // good enough
4839    }
4840  }
4841  int iColumn;
4842  const int * columnLength = matrix_->getVectorLengths();
4843  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4844    if (columnLowerWork_[iColumn]<columnUpperWork_[iColumn]) {
4845      int length = columnLength[iColumn];
4846      if (length>2) {
4847        maxLength = CoinMax(maxLength,length);
4848        minLength = CoinMin(minLength,length);
4849      }
4850    }
4851  }
4852  // If > 70 then do rows
4853  if (perturbation_>=70) {
4854    modifyRowCosts=true;
4855    perturbation_ -= 20;
4856    printf("Row costs modified, ");
4857  }
4858  bool uniformChange=false;
4859  bool inCbcOrOther = (specialOptions_&0x03000000)!=0;
4860  if (perturbation_>50) {
4861    // Experiment
4862    // maximumFraction could be 1.0e-10 to 1.0
4863    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};
4864    int whichOne = perturbation_-51;
4865    //if (inCbcOrOther&&whichOne>0)
4866    //whichOne--;
4867    maximumFraction = m[CoinMin(whichOne,10)];
4868  } else if (inCbcOrOther) {
4869    //maximumFraction = 1.0e-6;
4870  }
4871  int iRow;
4872  double smallestNonZero=1.0e100;
4873  numberNonZero=0;
4874  if (perturbation_>=50) {
4875    perturbation = 1.0e-8;
4876    bool allSame=true;
4877    double lastValue=0.0;
4878    for (iRow=0;iRow<numberRows_;iRow++) {
4879      double lo = rowLowerWork_[iRow];
4880      double up = rowUpperWork_[iRow];
4881      if (lo<up) {
4882        double value = fabs(rowObjectiveWork_[iRow]);
4883        perturbation = CoinMax(perturbation,value);
4884        if (value) {
4885          modifyRowCosts=true;
4886          smallestNonZero = CoinMin(smallestNonZero,value);
4887        }
4888      } 
4889      if (lo&&lo>-1.0e10) {
4890        numberNonZero++;
4891        lo=fabs(lo);
4892        if (!lastValue) 
4893          lastValue=lo;
4894        else if (fabs(lo-lastValue)>1.0e-7)
4895          allSame=false;
4896      }
4897      if (up&&up<1.0e10) {
4898        numberNonZero++;
4899        up=fabs(up);
4900        if (!lastValue) 
4901          lastValue=up;
4902        else if (fabs(up-lastValue)>1.0e-7)
4903          allSame=false;
4904      }
4905    }
4906    double lastValue2=0.0;
4907    for (iColumn=0;iColumn<numberColumns_;iColumn++) { 
4908      double lo = columnLowerWork_[iColumn];
4909      double up = columnUpperWork_[iColumn];
4910      if (lo<up) {
4911        double value = 
4912          fabs(objectiveWork_[iColumn]);
4913        perturbation = CoinMax(perturbation,value);
4914        if (value) {
4915          smallestNonZero = CoinMin(smallestNonZero,value);
4916        }
4917      }
4918      if (lo&&lo>-1.0e10) {
4919        //numberNonZero++;
4920        lo=fabs(lo);
4921        if (!lastValue2) 
4922          lastValue2=lo;
4923        else if (fabs(lo-lastValue2)>1.0e-7)
4924          allSame=false;
4925      }
4926      if (up&&up<1.0e10) {
4927        //numberNonZero++;
4928        up=fabs(up);
4929        if (!lastValue2) 
4930          lastValue2=up;
4931        else if (fabs(up-lastValue2)>1.0e-7)
4932          allSame=false;
4933      }
4934    }
4935    if (allSame) {
4936      // Check elements
4937      double smallestNegative;
4938      double largestNegative;
4939      double smallestPositive;
4940      double largestPositive;
4941      matrix_->rangeOfElements(smallestNegative,largestNegative,
4942                               smallestPositive,largestPositive);
4943      if (smallestNegative==largestNegative&&
4944          smallestPositive==largestPositive) {
4945        // Really hit perturbation
4946        double adjust = CoinMin(100.0*maximumFraction,1.0e-3*CoinMax(lastValue,lastValue2));
4947        maximumFraction = CoinMax(adjust,maximumFraction);
4948      }
4949    }
4950    perturbation = CoinMin(perturbation,smallestNonZero/maximumFraction);
4951  } else {
4952    // user is in charge
4953    maximumFraction = 1.0e-1;
4954    // but some experiments
4955    if (perturbation_<=-900) {
4956      modifyRowCosts=true;
4957      perturbation_ += 1000;
4958      printf("Row costs modified, ");
4959    }
4960    if (perturbation_<=-10) {
4961      perturbation_ += 10; 
4962      maximumFraction = 1.0;
4963      if ((-perturbation_)%100>=10) {
4964        uniformChange=true;
4965        perturbation_+=20;
4966      }
4967      while (perturbation_<-10) {
4968        perturbation_ += 100;
4969        maximumFraction *= 1.0e-1;
4970      }
4971    }
4972    perturbation = pow(10.0,perturbation_);
4973  }
4974  double largestZero=0.0;
4975  double largest=0.0;
4976  double largestPerCent=0.0;
4977  // modify costs
4978  bool printOut=(handler_->logLevel()==63);
4979  printOut=false;
4980  if (modifyRowCosts) {
4981    for (iRow=0;iRow<numberRows_;iRow++) {
4982      if (rowLowerWork_[iRow]<rowUpperWork_[iRow]) {
4983        double value = perturbation;
4984        double currentValue = rowObjectiveWork_[iRow];
4985        value = CoinMin(value,maximumFraction*(fabs(currentValue)+1.0e-1*perturbation+1.0e-3));
4986        if (rowLowerWork_[iRow]>-largeValue_) {
4987          if (fabs(rowLowerWork_[iRow])<fabs(rowUpperWork_[iRow])) 
4988            value *= randomNumberGenerator_.randomDouble();
4989          else
4990            value *= -randomNumberGenerator_.randomDouble();
4991        } else if (rowUpperWork_[iRow]<largeValue_) {
4992          value *= -randomNumberGenerator_.randomDouble();
4993        } else {
4994          value=0.0;
4995        }
4996        if (currentValue) {
4997          largest = CoinMax(largest,fabs(value));
4998          if (fabs(value)>fabs(currentValue)*largestPerCent)
4999            largestPerCent=fabs(value/currentValue);
5000        } else {
5001          largestZero=CoinMax(largestZero,fabs(value));
5002        }
5003        if (printOut)
5004          printf("row %d cost %g change %g\n",iRow,rowObjectiveWork_[iRow],value);
5005        rowObjectiveWork_[iRow] += value;
5006      }
5007    }
5008  }
5009  // 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};
5010  // 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};
5011  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};
5012  //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};
5013  double extraWeight=10.0;
5014  // Scale back if wanted
5015  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};
5016  if (constantPerturbation<99.0*dualTolerance_) {
5017    perturbation *= 0.1;
5018    extraWeight=0.5;
5019    memcpy(weight,weight2,sizeof(weight2));
5020  }
5021  // adjust weights if all columns long
5022  double factor=1.0;
5023  if (maxLength) {
5024    factor = 3.0/(double) minLength;
5025  }
5026  // Make variables with more elements more expensive
5027  const double m1 = 0.5;
5028  double smallestAllowed = CoinMin(1.0e-2*dualTolerance_,maximumFraction);
5029  //double largestAllowed = CoinMax(1.0e3*dualTolerance_,maximumFraction*10.0*averageCost);
5030  double largestAllowed = CoinMax(1.0e3*dualTolerance_,maximumFraction*averageCost);
5031  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
5032    if (columnLowerWork_[iColumn]<columnUpperWork_[iColumn]&&getStatus(iColumn)!=basic) {
5033      double value = perturbation;
5034      double currentValue = objectiveWork_[iColumn];
5035      value = CoinMin(value,constantPerturbation+maximumFraction*(fabs(currentValue)+1.0e-1*perturbation+1.0e-8));
5036      //value = CoinMin(value,constantPerturbation;+maximumFraction*fabs(currentValue));
5037      double value2 = constantPerturbation+1.0e-1*smallestNonZero;
5038      if (uniformChange) {
5039        value = maximumFraction;
5040        value2=maximumFraction;
5041      }
5042      if (columnLowerWork_[iColumn]>-largeValue_) {
5043        if (fabs(columnLowerWork_[iColumn])<
5044            fabs(columnUpperWork_[iColumn])) {
5045          value *= (1.0-m1 +m1*randomNumberGenerator_.randomDouble());
5046          value2 *= (1.0-m1 +m1*randomNumberGenerator_.randomDouble());
5047        } else {
5048          //value *= -(1.0-m1+m1*randomNumberGenerator_.randomDouble());
5049          //value2 *= -(1.0-m1+m1*randomNumberGenerator_.randomDouble());
5050          value=0.0;
5051        }
5052      } else if (columnUpperWork_[iColumn]<largeValue_) {
5053        value *= -(1.0-m1+m1*randomNumberGenerator_.randomDouble());
5054        value2 *= -(1.0-m1+m1*randomNumberGenerator_.randomDouble());
5055      } else {
5056        value=0.0;
5057      }
5058      if (value) {
5059        int length = columnLength[iColumn];
5060        if (length>3) {
5061          length = (int) ((double) length * factor);
5062          length = CoinMax(3,length);
5063        }
5064        double multiplier;
5065#if 1
5066        if (length<10)
5067          multiplier=weight[length];
5068        else
5069          multiplier = weight[10];
5070#else
5071        if (length<10)
5072          multiplier=weight[length];
5073        else
5074          multiplier = weight[10]+extraWeight*(length-10);
5075        multiplier *= 0.5;
5076#endif
5077        value *= multiplier;
5078        value = CoinMin(value,value2);
5079        if (savePerturbation<50||savePerturbation>60) {
5080          if (fabs(value)<=dualTolerance_)
5081            value=0.0;
5082        } else if (value) {
5083          // get in range
5084          if (fabs(value)<=smallestAllowed) {
5085            value *= 10.0;
5086            while (fabs(value)<=smallestAllowed) 
5087              value *= 10.0;
5088          } else if (fabs(value)>largestAllowed) {
5089            value *= 0.1;
5090            while (fabs(value)>largestAllowed) 
5091              value *= 0.1;
5092          }
5093        }
5094        if (currentValue) {
5095          largest = CoinMax(largest,fabs(value));
5096          if (fabs(value)>fabs(currentValue)*largestPerCent)
5097            largestPerCent=fabs(value/currentValue);
5098        } else {
5099          largestZero=CoinMax(largestZero,fabs(value));
5100        }
5101        // but negative if at ub
5102        if (getStatus(iColumn)==atUpperBound)
5103          value = -value;
5104        if (printOut)
5105          printf("col %d cost %g change %g\n",iColumn,objectiveWork_[iColumn],value);
5106        objectiveWork_[iColumn] += value;
5107      }
5108    }
5109  }
5110  handler_->message(CLP_SIMPLEX_PERTURB,messages_)
5111    <<100.0*maximumFraction<<perturbation<<largest<<100.0*largestPerCent<<largestZero
5112    <<CoinMessageEol;
5113  // and zero changes
5114  //int nTotal = numberRows_+numberColumns_;
5115  //CoinZeroN(cost_+nTotal,nTotal);
5116  // say perturbed
5117  perturbation_=101;
5118  return 0;
5119}
5120/* For strong branching.  On input lower and upper are new bounds
5121   while on output they are change in objective function values
5122   (>1.0e50 infeasible).
5123   Return code is 0 if nothing interesting, -1 if infeasible both
5124   ways and +1 if infeasible one way (check values to see which one(s))
5125   Returns -2 if bad factorization
5126*/
5127int ClpSimplexDual::strongBranching(int numberVariables,const int * variables,
5128                                    double * newLower, double * newUpper,
5129                                    double ** outputSolution,
5130                                    int * outputStatus, int * outputIterations,
5131                                    bool stopOnFirstInfeasible,
5132                                    bool alwaysFinish,
5133                                    int startFinishOptions)
5134{
5135  int i;
5136  int returnCode=0;
5137  double saveObjectiveValue = objectiveValue_;
5138  algorithm_ = -1;
5139
5140  //scaling(false);
5141
5142  // put in standard form (and make row copy)
5143  // create modifiable copies of model rim and do optional scaling
5144  createRim(7+8+16+32,true,startFinishOptions);
5145
5146  // change newLower and newUpper if scaled
5147
5148  // Do initial factorization
5149  // and set certain stuff
5150  // We can either set increasing rows so ...IsBasic gives pivot row
5151  // or we can just increment iBasic one by one
5152  // for now let ...iBasic give pivot row
5153  int useFactorization=false;
5154  if ((startFinishOptions&2)!=0&&(whatsChanged_&(2+512))==2+512)
5155    useFactorization=true; // Keep factorization if possible
5156  // switch off factorization if bad
5157  if (pivotVariable_[0]<0)
5158    useFactorization=false;
5159  if (!useFactorization||factorization_->numberRows()!=numberRows_) {
5160    useFactorization = false;
5161    factorization_->setDefaultValues();
5162
5163    int factorizationStatus = internalFactorize(0);
5164    if (factorizationStatus<0) {
5165      // some error
5166      // we should either debug or ignore
5167#ifndef NDEBUG
5168      printf("***** ClpDual strong branching factorization error - debug\n");
5169#endif
5170      return -2;
5171    } else if (factorizationStatus&&factorizationStatus<=numberRows_) {
5172      handler_->message(CLP_SINGULARITIES,messages_)
5173        <<factorizationStatus
5174        <<CoinMessageEol;
5175    }
5176  }
5177  // save stuff
5178  ClpFactorization saveFactorization(*factorization_);
5179  // save basis and solution
5180  double * saveSolution = new double[numberRows_+numberColumns_];
5181  CoinMemcpyN(solution_,
5182         numberRows_+numberColumns_,saveSolution);
5183  unsigned char * saveStatus = 
5184    new unsigned char [numberRows_+numberColumns_];
5185  CoinMemcpyN(status_,numberColumns_+numberRows_,saveStatus);
5186  // save bounds as createRim makes clean copies
5187  double * saveLower = new double[numberRows_+numberColumns_];
5188  CoinMemcpyN(lower_,
5189         numberRows_+numberColumns_,saveLower);
5190  double * saveUpper = new double[numberRows_+numberColumns_];
5191  CoinMemcpyN(upper_,
5192         numberRows_+numberColumns_,saveUpper);
5193  double * saveObjective = new double[numberRows_+numberColumns_];
5194  CoinMemcpyN(cost_,
5195         numberRows_+numberColumns_,saveObjective);
5196  int * savePivot = new int [numberRows_];
5197  CoinMemcpyN(pivotVariable_, numberRows_,savePivot);
5198  // need to save/restore weights.
5199
5200  int iSolution = 0;
5201  for (i=0;i<numberVariables;i++) {
5202    int iColumn = variables[i];
5203    double objectiveChange;
5204    double saveBound;
5205   
5206    // try down
5207   
5208    saveBound = columnUpper_[iColumn];
5209    // external view - in case really getting optimal
5210    columnUpper_[iColumn] = newUpper[i];
5211    if (scalingFlag_<=0) 
5212      upper_[iColumn] = newUpper[i]*rhsScale_;
5213    else 
5214      upper_[iColumn] = (newUpper[i]/columnScale_[iColumn])*rhsScale_; // scale
5215    // Start of fast iterations
5216    int status = fastDual(alwaysFinish);
5217    CoinAssert (problemStatus_||objectiveValue_<1.0e50);
5218    // make sure plausible
5219    double obj = CoinMax(objectiveValue_,saveObjectiveValue);
5220    if (status&&problemStatus_!=3) {
5221      // not finished - might be optimal
5222      checkPrimalSolution(rowActivityWork_,columnActivityWork_);
5223      double limit = 0.0;
5224      getDblParam(ClpDualObjectiveLimit, limit);
5225      if (!numberPrimalInfeasibilities_&&obj<limit) { 
5226        problemStatus_=0;
5227      } 
5228      status=problemStatus_;
5229    }
5230    if (status||(problemStatus_==0&&!isDualObjectiveLimitReached())) {
5231      objectiveChange = obj-saveObjectiveValue;
5232    } else {
5233      objectiveChange = 1.0e100;
5234      status=1;
5235    }
5236    if (problemStatus_==3)
5237      status=2;
5238
5239    if (scalingFlag_<=0) {
5240      CoinMemcpyN(solution_,numberColumns_,outputSolution[iSolution]);
5241    } else {
5242      int j;
5243      double * sol = outputSolution[iSolution];
5244      for (j=0;j<numberColumns_;j++) 
5245        sol[j] = solution_[j]*columnScale_[j];
5246    }
5247    outputStatus[iSolution]=status;
5248    outputIterations[iSolution]=numberIterations_;
5249    iSolution++;
5250    // restore
5251    CoinMemcpyN(saveSolution,
5252            numberRows_+numberColumns_,solution_);
5253    CoinMemcpyN(saveStatus,numberColumns_+numberRows_,status_);
5254    CoinMemcpyN(saveLower,
5255            numberRows_+numberColumns_,lower_);
5256    CoinMemcpyN(saveUpper,
5257            numberRows_+numberColumns_,upper_);
5258    CoinMemcpyN(saveObjective,
5259            numberRows_+numberColumns_,cost_);
5260    columnUpper_[iColumn] = saveBound;
5261    CoinMemcpyN(savePivot, numberRows_,pivotVariable_);
5262    delete factorization_;
5263    factorization_ = new ClpFactorization(saveFactorization,numberRows_);
5264
5265    newUpper[i]=objectiveChange;
5266#ifdef CLP_DEBUG
5267    printf("down on %d costs %g\n",iColumn,objectiveChange);
5268#endif
5269         
5270    // try up
5271   
5272    saveBound = columnLower_[iColumn];
5273    // external view - in case really getting optimal
5274    columnLower_[iColumn] = newLower[i];
5275    if (scalingFlag_<=0) 
5276      lower_[iColumn] = newLower[i]*rhsScale_;
5277    else 
5278      lower_[iColumn] = (newLower[i]/columnScale_[iColumn])*rhsScale_; // scale
5279    // Start of fast iterations
5280    status = fastDual(alwaysFinish);
5281    // make sure plausible
5282    obj = CoinMax(objectiveValue_,saveObjectiveValue);
5283    if (status&&problemStatus_!=3) {
5284      // not finished - might be optimal
5285      checkPrimalSolution(rowActivityWork_,columnActivityWork_);
5286      double limit = 0.0;
5287      getDblParam(ClpDualObjectiveLimit, limit);
5288      if (!numberPrimalInfeasibilities_&&obj< limit) { 
5289        problemStatus_=0;
5290      } 
5291      status=problemStatus_;
5292    }
5293    if (status||(problemStatus_==0&&!isDualObjectiveLimitReached())) {
5294      objectiveChange = obj-saveObjectiveValue;
5295    } else {
5296      objectiveChange = 1.0e100;
5297      status=1;
5298    }
5299    if (problemStatus_==3)
5300      status=2;
5301    if (scalingFlag_<=0) {
5302      CoinMemcpyN(solution_,numberColumns_,outputSolution[iSolution]);
5303    } else {
5304      int j;
5305      double * sol = outputSolution[iSolution];
5306      for (j=0;j<numberColumns_;j++) 
5307        sol[j] = solution_[j]*columnScale_[j];
5308    }
5309    outputStatus[iSolution]=status;
5310    outputIterations[iSolution]=numberIterations_;
5311    iSolution++;
5312
5313    // restore
5314    CoinMemcpyN(saveSolution,
5315            numberRows_+numberColumns_,solution_);
5316    CoinMemcpyN(saveStatus,numberColumns_+numberRows_,status_);
5317    CoinMemcpyN(saveLower,
5318            numberRows_+numberColumns_,lower_);
5319    CoinMemcpyN(saveUpper,
5320            numberRows_+numberColumns_,upper_);
5321    CoinMemcpyN(saveObjective,
5322            numberRows_+numberColumns_,cost_);
5323    columnLower_[iColumn] = saveBound;
5324    CoinMemcpyN(savePivot, numberRows_,pivotVariable_);
5325    delete factorization_;
5326    factorization_ = new ClpFactorization(saveFactorization,numberRows_);
5327
5328    newLower[i]=objectiveChange;
5329#ifdef CLP_DEBUG
5330    printf("up on %d costs %g\n",iColumn,objectiveChange);
5331#endif
5332         
5333    /* Possibilities are:
5334       Both sides feasible - store
5335       Neither side feasible - set objective high and exit if desired
5336       One side feasible - change bounds and resolve
5337    */
5338    if (newUpper[i]<1.0e100) {
5339      if(newLower[i]<1.0e100) {
5340        // feasible - no action
5341      } else {
5342        // up feasible, down infeasible
5343        returnCode=1;
5344        if (stopOnFirstInfeasible)
5345          break;
5346      }
5347    } else {
5348      if(newLower[i]<1.0e100) {
5349        // down feasible, up infeasible
5350        returnCode=1;
5351        if (stopOnFirstInfeasible)
5352          break;
5353      } else {
5354        // neither side feasible
5355        returnCode=-1;
5356        break;
5357      }
5358    }
5359  }
5360  delete [] saveSolution;
5361  delete [] saveLower;
5362  delete [] saveUpper;
5363  delete [] saveObjective;
5364  delete [] saveStatus;
5365  delete [] savePivot;
5366  if ((startFinishOptions&1)==0) {
5367    deleteRim(1);
5368    whatsChanged_=0;
5369  } else {
5370    // Original factorization will have been put back by last loop
5371    //delete factorization_;
5372    //factorization_ = new ClpFactorization(saveFactorization);
5373    deleteRim(0);
5374    // mark all as current
5375    whatsChanged_ = 0xffff;
5376  }
5377  objectiveValue_ = saveObjectiveValue;
5378  return returnCode;
5379}
5380// treat no pivot as finished (unless interesting)
5381int ClpSimplexDual::fastDual(bool alwaysFinish)
5382{
5383  algorithm_ = -1;
5384  secondaryStatus_=0;
5385  // Say in fast dual
5386  specialOptions_ |= 16384;
5387  int saveDont = dontFactorizePivots_;
5388  if ((specialOptions_&2048)==0)
5389    dontFactorizePivots_=0;
5390  else if(!dontFactorizePivots_)
5391    dontFactorizePivots_ = 20;
5392  //handler_->setLogLevel(63);
5393  // save data
5394  ClpDataSave data = saveData();
5395  dualTolerance_=dblParam_[ClpDualTolerance];
5396  primalTolerance_=dblParam_[ClpPrimalTolerance];
5397
5398  // save dual bound
5399  double saveDualBound = dualBound_;
5400
5401  if (alphaAccuracy_!=-1.0)
5402    alphaAccuracy_ = 1.0;
5403  double objectiveChange;
5404  // for dual we will change bounds using dualBound_
5405  // for this we need clean basis so it is after factorize
5406#if 0
5407  {
5408    int numberTotal = numberRows_+numberColumns_;
5409    double * saveSol = CoinCopyOfArray(solution_,numberTotal);
5410    double * saveDj = CoinCopyOfArray(dj_,numberTotal);
5411    double tolerance = 1.0e-8;
5412    gutsOfSolution(NULL,NULL);
5413    int j;
5414    double largestPrimal = tolerance;
5415    int iPrimal=-1;
5416    for (j=0;j<numberTotal;j++) {
5417      double difference = solution_[j]-saveSol[j];
5418      if (fabs(difference)>largestPrimal) {
5419        iPrimal=j;
5420        largestPrimal=fabs(difference);
5421      }
5422    }
5423    double largestDual = tolerance;
5424    int iDual=-1;
5425    for (j=0;j<numberTotal;j++) {
5426      double difference = dj_[j]-saveDj[j];
5427      if (fabs(difference)>largestDual&&upper_[j]>lower_[j]) {
5428        iDual=j;
5429        largestDual=fabs(difference);
5430      }
5431    }
5432    if (iPrimal>=0||iDual>=0)
5433        printf("pivots %d primal diff(%g,%d) dual diff(%g,%d)\n",
5434               factorization_->pivots(),
5435               largestPrimal,iPrimal,
5436               largestDual,iDual);
5437    delete [] saveSol;
5438    delete [] saveDj;
5439  }
5440#else
5441  if ((specialOptions_&524288)==0)
5442    gutsOfSolution(NULL,NULL);
5443#endif
5444#if 0
5445  if (numberPrimalInfeasibilities_!=1||
5446      numberDualInfeasibilities_)
5447    printf("dual %g (%d) primal %g (%d)\n",
5448           sumDualInfeasibilities_,numberDualInfeasibilities_,
5449           sumPrimalInfeasibilities_,numberPrimalInfeasibilities_);
5450#endif
5451  numberFake_ =0; // Number of variables at fake bounds
5452  numberChanged_ =0; // Number of variables with changed costs
5453  changeBounds(1,NULL,objectiveChange);
5454
5455  problemStatus_ = -1;
5456  numberIterations_=0;
5457  if ((specialOptions_&524288)==0) {
5458    factorization_->sparseThreshold(0);
5459    factorization_->goSparse();
5460  }
5461
5462  int lastCleaned=0; // last time objective or bounds cleaned up
5463
5464  // number of times we have declared optimality
5465  numberTimesOptimal_=0;
5466
5467  // This says whether to restore things etc
5468  int factorType=0;
5469  /*
5470    Status of problem:
5471    0 - optimal
5472    1 - infeasible
5473    2 - unbounded
5474    -1 - iterating
5475    -2 - factorization wanted
5476    -3 - redo checking without factorization
5477    -4 - looks infeasible
5478
5479    BUT also from whileIterating return code is:
5480
5481   -1 iterations etc
5482   -2 inaccuracy
5483   -3 slight inaccuracy (and done iterations)
5484   +0 looks optimal (might be unbounded - but we will investigate)
5485   +1 looks infeasible
5486   +3 max iterations
5487
5488  */
5489
5490  int returnCode = 0;
5491
5492  int iRow,iColumn;
5493  int maxPass = maximumIterations();
5494  while (problemStatus_<0) {
5495    // clear
5496    for (iRow=0;iRow<4;iRow++) {
5497      rowArray_[iRow]->clear();
5498    }   
5499   
5500    for (iColumn=0;iColumn<2;iColumn++) {
5501      columnArray_[iColumn]->clear();
5502    }   
5503
5504    // give matrix (and model costs and bounds a chance to be
5505    // refreshed (normally null)
5506    matrix_->refresh(this);
5507    // If getting nowhere - why not give it a kick
5508    // does not seem to work too well - do some more work
5509    if ((specialOptions_&524288)!=0&&
5510        perturbation_<101&&numberIterations_>2*(numberRows_+numberColumns_)) {
5511      perturb();
5512      // Can't get here if values pass
5513      gutsOfSolution(NULL,NULL);
5514      if (handler_->logLevel()>2) {
5515        handler_->message(CLP_SIMPLEX_STATUS,messages_)
5516          <<numberIterations_<<objectiveValue();
5517        handler_->printing(sumPrimalInfeasibilities_>0.0)
5518          <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
5519        handler_->printing(sumDualInfeasibilities_>0.0)
5520          <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
5521        handler_->printing(numberDualInfeasibilitiesWithoutFree_
5522                           <numberDualInfeasibilities_)
5523                             <<numberDualInfeasibilitiesWithoutFree_;
5524        handler_->message()<<CoinMessageEol;
5525      }
5526    }
5527    // may factorize, checks if problem finished
5528    // should be able to speed this up on first time
5529    statusOfProblemInDual(lastCleaned,factorType,NULL,data,0);
5530
5531    // Say good factorization
5532    factorType=1;
5533    maxPass--;
5534    if (maxPass<-10) {
5535      // odd
5536      returnCode=1;
5537      problemStatus_ = 3; 
5538      // can't say anything interesting - might as well return
5539#ifdef CLP_DEBUG
5540      printf("returning from fastDual after %d iterations with code %d because of loop\n",
5541             numberIterations_,returnCode);
5542#endif
5543      break;
5544    }
5545
5546    // Do iterations
5547    if (problemStatus_<0) {
5548      double * givenPi=NULL;
5549      returnCode = whileIterating(givenPi,0);
5550      if ((!alwaysFinish&&returnCode<0)||returnCode==3) {
5551        if (returnCode!=3)
5552          assert (problemStatus_<0);
5553        returnCode=1;
5554        problemStatus_ = 3; 
5555        // can't say anything interesting - might as well return
5556#ifdef CLP_DEBUG
5557        printf("returning from fastDual after %d iterations with code %d\n",
5558               numberIterations_,returnCode);
5559#endif
5560        break;
5561      }
5562      if (returnCode==-2)
5563        factorType=3;
5564      returnCode=0;
5565    }
5566  }
5567
5568  // clear
5569  for (iRow=0;iRow<4;iRow++) {
5570    rowArray_[iRow]->clear();
5571  }   
5572 
5573  for (iColumn=0;iColumn<2;iColumn++) {
5574    columnArray_[iColumn]->clear();
5575  }   
5576  // Say not in fast dual
5577  specialOptions_ &= ~16384;
5578  assert(!numberFake_||((specialOptions_&(2048|4096))!=0&&dualBound_>=1.0e8)
5579         ||returnCode||problemStatus_); // all bounds should be okay
5580  if (numberFake_>0) {
5581    // Set back
5582    double dummy;
5583    changeBounds(2,NULL,dummy);
5584  }
5585  // Restore any saved stuff
5586  restoreData(data);
5587  dontFactorizePivots_ = saveDont;
5588  dualBound_ = saveDualBound;
5589  return returnCode;
5590}
5591// This does first part of StrongBranching
5592ClpFactorization * 
5593ClpSimplexDual::setupForStrongBranching(char * arrays, int numberRows, int numberColumns)
5594{
5595  algorithm_ = -1;
5596  // put in standard form (and make row copy)
5597  // create modifiable copies of model rim and do optional scaling
5598  int startFinishOptions;
5599  /*  COIN_CLP_VETTED
5600      Looks safe for Cbc
5601  */
5602  bool safeOption = ((specialOptions_&COIN_CBC_USING_CLP)!=0);
5603  if((specialOptions_&4096)==0||(auxiliaryModel_&&!safeOption)) {
5604    startFinishOptions=0;
5605  } else {
5606    startFinishOptions=1+2+4;
5607  }
5608  createRim(7+8+16+32,true,startFinishOptions);
5609  // Do initial factorization
5610  // and set certain stuff
5611  // We can either set increasing rows so ...IsBasic gives pivot row
5612  // or we can just increment iBasic one by one
5613  // for now let ...iBasic give pivot row
5614  bool useFactorization=false;
5615  if ((startFinishOptions&2)!=0&&(whatsChanged_&(2+512))==2+512) {
5616    useFactorization=true; // Keep factorization if possible
5617    // switch off factorization if bad
5618    if (pivotVariable_[0]<0||factorization_->numberRows()!=numberRows_) 
5619      useFactorization=false;
5620  }
5621  if (!useFactorization) {
5622    factorization_->setDefaultValues();
5623
5624    int factorizationStatus = internalFactorize(0);
5625    if (factorizationStatus<0) {
5626      // some error
5627      // we should either debug or ignore
5628#ifndef NDEBUG
5629      printf("***** ClpDual strong branching factorization error - debug\n");
5630#endif
5631    } else if (factorizationStatus&&factorizationStatus<=numberRows_) {
5632      handler_->message(CLP_SINGULARITIES,messages_)
5633        <<factorizationStatus
5634        <<CoinMessageEol;
5635    }
5636  }
5637  double * arrayD = (double *) arrays;
5638  arrayD[0]=objectiveValue()*optimizationDirection_;
5639  double * saveSolution = arrayD+1;
5640  double * saveLower = saveSolution + (numberRows+numberColumns);
5641  double * saveUpper = saveLower + (numberRows+numberColumns);
5642  double * saveObjective = saveUpper + (numberRows+numberColumns);
5643  double * saveLowerOriginal = saveObjective + (numberRows+numberColumns);
5644  double * saveUpperOriginal = saveLowerOriginal + numberColumns;
5645  arrayD = saveUpperOriginal + numberColumns;
5646  int * savePivot = (int *) arrayD;
5647  int * whichRow = savePivot+numberRows;
5648  int * whichColumn = whichRow + 3*numberRows;
5649  int * arrayI = whichColumn + 2*numberColumns;
5650  unsigned char * saveStatus = (unsigned char *) (arrayI+1);
5651  // save stuff
5652  // save basis and solution
5653  CoinMemcpyN(solution_,
5654          numberRows_+numberColumns_,saveSolution);
5655  CoinMemcpyN(status_,numberColumns_+numberRows_,saveStatus);
5656  CoinMemcpyN(lower_,
5657          numberRows_+numberColumns_,saveLower);
5658  CoinMemcpyN(upper_,
5659          numberRows_+numberColumns_,saveUpper);
5660  CoinMemcpyN(cost_,
5661          numberRows_+numberColumns_,saveObjective);
5662  CoinMemcpyN(pivotVariable_, numberRows_,savePivot);
5663  return new ClpFactorization(*factorization_,numberRows_);
5664}
5665// This cleans up after strong branching
5666void 
5667ClpSimplexDual::cleanupAfterStrongBranching()
5668{
5669  int startFinishOptions;
5670  /*  COIN_CLP_VETTED
5671      Looks safe for Cbc
5672  */
5673  bool safeOption = ((specialOptions_&COIN_CBC_USING_CLP)!=0);
5674  if((specialOptions_&4096)==0||(auxiliaryModel_&&!safeOption)) {
5675    startFinishOptions=0;
5676  } else {
5677    startFinishOptions=1+2+4;
5678  }
5679  if ((startFinishOptions&1)==0) {
5680    deleteRim(1);
5681    whatsChanged_=0;
5682  } else {
5683    // Original factorization will have been put back by last loop
5684    //delete factorization_;
5685    //factorization_ = new ClpFactorization(saveFactorization);
5686    deleteRim(0);
5687    // mark all as current
5688    whatsChanged_ = 0xffff;
5689  }
5690}
5691/* Checks number of variables at fake bounds.  This is used by fastDual
5692   so can exit gracefully before end */
5693int 
5694ClpSimplexDual::numberAtFakeBound()
5695{
5696  int iSequence;
5697  int numberFake=0;
5698 
5699  for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
5700    FakeBound bound = getFakeBound(iSequence);
5701    switch(getStatus(iSequence)) {
5702
5703    case basic:
5704      break;
5705    case isFree:
5706    case superBasic:
5707    case ClpSimplex::isFixed:
5708      //setFakeBound (iSequence, noFake);
5709      break;
5710    case atUpperBound:
5711      if (bound==upperFake||bound==bothFake)
5712        numberFake++;
5713      break;
5714    case atLowerBound:
5715      if (bound==lowerFake||bound==bothFake)
5716        numberFake++;
5717      break;
5718    }
5719  }
5720  numberFake_ = numberFake;
5721  return numberFake;
5722}
5723/* Pivot out a variable and choose an incoing one.  Assumes dual
5724   feasible - will not go through a reduced cost. 
5725   Returns step length in theta
5726   Returns ray in ray_ (or NULL if no pivot)
5727   Return codes as before but -1 means no acceptable pivot
5728*/
5729int 
5730ClpSimplexDual::pivotResult()
5731{
5732  abort();
5733  return 0;
5734}
5735/*
5736   Row array has row part of pivot row
5737   Column array has column part.
5738   This is used in dual values pass
5739*/
5740void
5741ClpSimplexDual::checkPossibleValuesMove(CoinIndexedVector * rowArray,
5742                                        CoinIndexedVector * columnArray,
5743                                        double acceptablePivot)
5744{
5745  double * work;
5746  int number;
5747  int * which;
5748  int iSection;
5749
5750  double tolerance = dualTolerance_*1.001;
5751
5752  double thetaDown = 1.0e31;
5753  double changeDown ;
5754  double thetaUp = 1.0e31;
5755  double bestAlphaDown = acceptablePivot*0.99999;
5756  double bestAlphaUp = acceptablePivot*0.99999;
5757  int sequenceDown =-1;
5758  int sequenceUp = sequenceOut_;
5759
5760  double djBasic = dj_[sequenceOut_];
5761  if (djBasic>0.0) {
5762    // basic at lower bound so directionOut_ 1 and -1 in pivot row
5763    // dj will go to zero on other way
5764    thetaUp = djBasic;
5765    changeDown = -lower_[sequenceOut_];
5766  } else {
5767    // basic at upper bound so directionOut_ -1 and 1 in pivot row
5768    // dj will go to zero on other way
5769    thetaUp = -djBasic;
5770    changeDown = upper_[sequenceOut_];
5771  }
5772  bestAlphaUp = 1.0;
5773  int addSequence;
5774
5775  double alphaUp=0.0;
5776  double alphaDown=0.0;
5777
5778  for (iSection=0;iSection<2;iSection++) {
5779
5780    int i;
5781    if (!iSection) {
5782      work = rowArray->denseVector();
5783      number = rowArray->getNumElements();
5784      which = rowArray->getIndices();
5785      addSequence = numberColumns_;
5786    } else {
5787      work = columnArray->denseVector();
5788      number = columnArray->getNumElements();
5789      which = columnArray->getIndices();
5790      addSequence = 0;
5791    }
5792   
5793    for (i=0;i<number;i++) {
5794      int iSequence = which[i];
5795      int iSequence2 = iSequence + addSequence;
5796      double alpha;
5797      double oldValue;
5798      double value;
5799
5800      switch(getStatus(iSequence2)) {
5801         
5802      case basic: 
5803        break;
5804      case ClpSimplex::isFixed:
5805        alpha = work[i];
5806        changeDown += alpha*upper_[iSequence2];
5807        break;
5808      case isFree:
5809      case superBasic:
5810        alpha = work[i];
5811        // dj must be effectively zero as dual feasible
5812        if (fabs(alpha)>bestAlphaUp) {
5813          thetaDown = 0.0;
5814          thetaUp = 0.0;
5815          bestAlphaDown = fabs(alpha);
5816          bestAlphaUp = bestAlphaDown;
5817          sequenceDown =iSequence2;
5818          sequenceUp = sequenceDown;
5819          alphaUp = alpha;
5820          alphaDown = alpha;
5821        }
5822        break;
5823      case atUpperBound:
5824        alpha = work[i];
5825        oldValue = dj_[iSequence2];
5826        changeDown += alpha*upper_[iSequence2];
5827        if (alpha>=acceptablePivot) {
5828          // might do other way
5829          value = oldValue+thetaUp*alpha;
5830          if (value>-tolerance) {
5831            if (value>tolerance||fabs(alpha)>bestAlphaUp) {
5832              thetaUp = -oldValue/alpha;
5833              bestAlphaUp = fabs(alpha);
5834              sequenceUp = iSequence2;
5835              alphaUp = alpha;
5836            }
5837          }
5838        } else if (alpha<=-acceptablePivot) {
5839          // might do this way
5840          value = oldValue-thetaDown*alpha;
5841          if (value>-tolerance) {
5842            if (value>tolerance||fabs(alpha)>bestAlphaDown) {
5843              thetaDown = oldValue/alpha;
5844              bestAlphaDown = fabs(alpha);
5845              sequenceDown = iSequence2;
5846              alphaDown = alpha;
5847            }
5848          }
5849        }
5850        break;
5851      case atLowerBound:
5852        alpha = work[i];
5853        oldValue = dj_[iSequence2];
5854        changeDown += alpha*lower_[iSequence2];
5855        if (alpha<=-acceptablePivot) {
5856          // might do other way
5857          value = oldValue+thetaUp*alpha;
5858          if (value<tolerance) {
5859            if (value<-tolerance||fabs(alpha)>bestAlphaUp) {
5860              thetaUp = -oldValue/alpha;
5861              bestAlphaUp = fabs(alpha);
5862              sequenceUp = iSequence2;
5863              alphaUp = alpha;
5864            }
5865          }
5866        } else if (alpha>=acceptablePivot) {
5867          // might do this way
5868          value = oldValue-thetaDown*alpha;
5869          if (value<tolerance) {
5870            if (value<-tolerance||fabs(alpha)>bestAlphaDown) {
5871              thetaDown = oldValue/alpha;
5872              bestAlphaDown = fabs(alpha);
5873              sequenceDown = iSequence2;
5874              alphaDown = alpha;
5875            }
5876          }
5877        }
5878        break;
5879      }
5880    }
5881  }