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

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

get rid of compiler warnings

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