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

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

ifdef some printing

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