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

Last change on this file since 1111 was 1111, checked in by forrest, 12 years ago

trying to improve networks

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