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

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

memory leak

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