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

Last change on this file since 1426 was 1426, checked in by forrest, 10 years ago

minor changes for an idea

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