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

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

changes for factorization and aux region

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