source: trunk/ClpSimplexDual.cpp @ 111

Last change on this file since 111 was 111, checked in by forrest, 17 years ago

Changes to make OsiSimplex? stuff easier

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 82.0 KB
Line 
1
2// Copyright (C) 2002, International Business Machines
3// Corporation and others.  All Rights Reserved.
4
5
6/* Notes on implementation of dual simplex algorithm.
7
8   When dual feasible:
9
10   If primal feasible, we are optimal.  Otherwise choose an infeasible
11   basic variable to leave basis (normally going to nearest bound) (B).  We
12   now need to find an incoming variable which will leave problem
13   dual feasible so we get the row of the tableau corresponding to
14   the basic variable (with the correct sign depending if basic variable
15   above or below feasibility region - as that affects whether reduced
16   cost on outgoing variable has to be positive or negative).
17
18   We now perform a ratio test to determine which incoming variable will
19   preserve dual feasibility (C).  If no variable found then problem
20   is infeasible (in primal sense).  If there is a variable, we then
21   perform pivot and repeat.  Trivial?
22
23   -------------------------------------------
24
25   A) How do we get dual feasible?  If all variables have bounds then
26   it is trivial to get feasible by putting non-basic variables to
27   correct bounds.  OSL did not have a phase 1/phase 2 approach but
28   instead effectively put fake bounds on variables and this is the
29   approach here, although I had hoped to make it cleaner.
30
31   If there is a weight of X on getting dual feasible:
32     Non-basic variables with negative reduced costs are put to
33     lesser of their upper bound and their lower bound + X.
34     Similarly, mutatis mutandis, for positive reduced costs.
35
36   Free variables should normally be in basis, otherwise I have
37   coding which may be able to come out (and may not be correct).
38
39   In OSL, this weight was changed heuristically, here at present
40   it is only increased if problem looks finished.  If problem is
41   feasible I check for unboundedness.  If not unbounded we
42   could play with going into primal.  As long as weights increase
43   any algorithm would be finite.
44   
45   B) Which outgoing variable to choose is a virtual base class.
46   For difficult problems steepest edge is preferred while for
47   very easy (large) problems we will need partial scan.
48
49   C) Sounds easy, but this is hardest part of algorithm.
50      1) Instead of stopping at first choice, we may be able
51      to flip that variable to other bound and if objective
52      still improving choose again.  These mini iterations can
53      increase speed by orders of magnitude but we may need to
54      go to more of a bucket choice of variable rather than looking
55      at them one by one (for speed).
56      2) Accuracy.  Reduced costs may be of wrong sign but less than
57      tolerance.  Pivoting on these makes objective go backwards.
58      OSL modified cost so a zero move was made, Gill et al
59      (in primal analogue) modified so a strictly positive move was
60      made.  It is not quite as neat in dual but that is what we
61      try and do.  The two problems are that re-factorizations can
62      change reduced costs above and below tolerances and that when
63      finished we need to reset costs and try again.
64      3) Degeneracy.  Gill et al helps but may not be enough.  We
65      may need more.  Also it can improve speed a lot if we perturb
66      the costs significantly. 
67
68  References:
69     Forrest and Goldfarb, Steepest-edge simplex algorithms for
70       linear programming - Mathematical Programming 1992
71     Forrest and Tomlin, Implementing the simplex method for
72       the Optimization Subroutine Library - IBM Systems Journal 1992
73     Gill, Murray, Saunders, Wright A Practical Anti-Cycling
74       Procedure for Linear and Nonlinear Programming SOL report 1988
75
76
77  TODO:
78 
79  a) Better recovery procedures.  At present I never check on forward
80     progress.  There is checkpoint/restart with reducing
81     re-factorization frequency, but this is only on singular
82     factorizations.
83  b) Fast methods for large easy problems (and also the option for
84     the code to automatically choose which method).
85  c) We need to be able to stop in various ways for OSI - this
86     is fairly easy.
87
88 */
89
90
91#include "CoinPragma.hpp"
92
93#include <math.h>
94
95#include "CoinHelperFunctions.hpp"
96#include "ClpSimplexDual.hpp"
97#include "ClpFactorization.hpp"
98#include "CoinPackedMatrix.hpp"
99#include "CoinIndexedVector.hpp"
100#include "CoinWarmStartBasis.hpp"
101#include "ClpDualRowDantzig.hpp"
102#include "ClpMessage.hpp"
103#include <cfloat>
104#include <cassert>
105#include <string>
106#include <stdio.h>
107#include <iostream>
108#define CHECK_DJ
109// dual
110int ClpSimplexDual::dual ( )
111{
112
113  /* *** Method
114     This is a vanilla version of dual simplex.
115
116     It tries to be a single phase approach with a weight of 1.0 being
117     given to getting optimal and a weight of dualBound_ being
118     given to getting dual feasible.  In this version I have used the
119     idea that this weight can be thought of as a fake bound.  If the
120     distance between the lower and upper bounds on a variable is less
121     than the feasibility weight then we are always better off flipping
122     to other bound to make dual feasible.  If the distance is greater
123     then we make up a fake bound dualBound_ away from one bound.
124     If we end up optimal or primal infeasible, we check to see if
125     bounds okay.  If so we have finished, if not we increase dualBound_
126     and continue (after checking if unbounded). I am undecided about
127     free variables - there is coding but I am not sure about it.  At
128     present I put them in basis anyway.
129
130     The code is designed to take advantage of sparsity so arrays are
131     seldom zeroed out from scratch or gone over in their entirety.
132     The only exception is a full scan to find outgoing variable.  This
133     will be changed to keep an updated list of infeasibilities (or squares
134     if steepest edge).  Also on easy problems we don't need full scan - just
135     pick first reasonable.
136
137     One problem is how to tackle degeneracy and accuracy.  At present
138     I am using the modification of costs which I put in OSL and which was
139     extended by Gill et al.  I am still not sure of the exact details.
140
141     The flow of dual is three while loops as follows:
142
143     while (not finished) {
144
145       while (not clean solution) {
146
147          Factorize and/or clean up solution by flipping variables so
148          dual feasible.  If looks finished check fake dual bounds.
149          Repeat until status is iterating (-1) or finished (0,1,2)
150
151       }
152
153       while (status==-1) {
154
155         Iterate until no pivot in or out or time to re-factorize.
156
157         Flow is:
158
159         choose pivot row (outgoing variable).  if none then
160         we are primal feasible so looks as if done but we need to
161         break and check bounds etc.
162
163         Get pivot row in tableau
164
165         Choose incoming column.  If we don't find one then we look
166         primal infeasible so break and check bounds etc.  (Also the
167         pivot tolerance is larger after any iterations so that may be
168         reason)
169
170         If we do find incoming column, we may have to adjust costs to
171         keep going forwards (anti-degeneracy).  Check pivot will be stable
172         and if unstable throw away iteration (we will need to implement
173         flagging of basic variables sometime) and break to re-factorize.
174         If minor error re-factorize after iteration.
175
176         Update everything (this may involve flipping variables to stay
177         dual feasible.
178
179       }
180
181     }
182
183     At present we never check we are going forwards.  I overdid that in
184     OSL so will try and make a last resort.
185
186     Needs partial scan pivot out option.
187     Needs dantzig, uninitialized and full steepest edge options (can still
188     use partial scan)
189
190     May need other anti-degeneracy measures, especially if we try and use
191     loose tolerances as a way to solve in fewer iterations.
192
193     I like idea of dynamic scaling.  This gives opportunity to decouple
194     different implications of scaling for accuracy, iteration count and
195     feasibility tolerance.
196
197  */
198
199  algorithm_ = -1;
200
201  // save data
202  ClpDataSave data = saveData();
203
204  // sanity check
205  // initialize - no values pass and algorithm_ is -1
206  // put in standard form (and make row copy)
207  // create modifiable copies of model rim and do optional scaling
208  // If problem looks okay
209  // Do initial factorization
210  // If user asked for perturbation - do it
211  if (!startup(0)) {
212    // looks okay
213   
214    double objectiveChange;
215    numberFake_ =0; // Number of variables at fake bounds
216    changeBounds(true,NULL,objectiveChange);
217   
218    int lastCleaned=0; // last time objective or bounds cleaned up
219   
220
221    // Progress indicator
222    ClpSimplexProgress progress(this);
223   
224    // This says whether to restore things etc
225    int factorType=0;
226    /*
227      Status of problem:
228      0 - optimal
229      1 - infeasible
230      2 - unbounded
231      -1 - iterating
232      -2 - factorization wanted
233      -3 - redo checking without factorization
234      -4 - looks infeasible
235    */
236    while (problemStatus_<0) {
237      int iRow, iColumn;
238      // clear
239      for (iRow=0;iRow<4;iRow++) {
240        rowArray_[iRow]->clear();
241      }   
242     
243      for (iColumn=0;iColumn<2;iColumn++) {
244        columnArray_[iColumn]->clear();
245      }   
246     
247      // give matrix (and model costs and bounds a chance to be
248      // refreshed (normally null)
249      matrix_->refresh(this);
250      // If getting nowhere - why not give it a kick
251#if 0
252      // does not seem to work too well - do some more work
253      if (perturbation_<101&&numberIterations_>2*(numberRows_+numberColumns_))
254        perturb();
255#endif
256      // may factorize, checks if problem finished
257      statusOfProblemInDual(lastCleaned,factorType,progress);
258     
259      // Say good factorization
260      factorType=1;
261      if (data.sparseThreshold_) {
262        // use default at present
263        factorization_->sparseThreshold(0);
264        factorization_->goSparse();
265      }
266     
267      // Do iterations
268      whileIterating();
269    }
270  }
271
272  assert (problemStatus_||!sumPrimalInfeasibilities_);
273
274  // clean up
275  finish();
276
277  // Restore any saved stuff
278  restoreData(data);
279  return problemStatus_;
280}
281/* Reasons to come out:
282   -1 iterations etc
283   -2 inaccuracy
284   -3 slight inaccuracy (and done iterations)
285   +0 looks optimal (might be unbounded - but we will investigate)
286   +1 looks infeasible
287   +3 max iterations
288 */
289int
290ClpSimplexDual::whileIterating()
291{
292#ifdef CLP_DEBUG
293  int debugIteration=-1;
294#endif
295  // if can't trust much and long way from optimal then relax
296  if (largestPrimalError_>10.0)
297    factorization_->relaxAccuracyCheck(min(1.0e2,largestPrimalError_/10.0));
298  else
299    factorization_->relaxAccuracyCheck(1.0);
300  // status stays at -1 while iterating, >=0 finished, -2 to invert
301  // status -3 to go to top without an invert
302  int returnCode = -1;
303  while (problemStatus_==-1) {
304#ifdef CLP_DEBUG
305    {
306      int i;
307      for (i=0;i<4;i++) {
308        rowArray_[i]->checkClear();
309      }   
310      for (i=0;i<2;i++) {
311        columnArray_[i]->checkClear();
312      }   
313    }     
314#endif
315#if CLP_DEBUG>2
316    // very expensive
317    if (numberIterations_>0&&numberIterations_<-801) {
318      handler_->setLogLevel(63);
319      double saveValue = objectiveValue_;
320      double * saveRow1 = new double[numberRows_];
321      double * saveRow2 = new double[numberRows_];
322      memcpy(saveRow1,rowReducedCost_,numberRows_*sizeof(double));
323      memcpy(saveRow2,rowActivityWork_,numberRows_*sizeof(double));
324      double * saveColumn1 = new double[numberColumns_];
325      double * saveColumn2 = new double[numberColumns_];
326      memcpy(saveColumn1,reducedCostWork_,numberColumns_*sizeof(double));
327      memcpy(saveColumn2,columnActivityWork_,numberColumns_*sizeof(double));
328      gutsOfSolution(rowActivityWork_,columnActivityWork_);
329      printf("xxx %d old obj %g, recomputed %g, sum dual inf %g\n",
330             numberIterations_,
331             saveValue,objectiveValue_,sumDualInfeasibilities_);
332      if (saveValue>objectiveValue_+1.0e-2)
333        printf("**bad**\n");
334      memcpy(rowReducedCost_,saveRow1,numberRows_*sizeof(double));
335      memcpy(rowActivityWork_,saveRow2,numberRows_*sizeof(double));
336      memcpy(reducedCostWork_,saveColumn1,numberColumns_*sizeof(double));
337      memcpy(columnActivityWork_,saveColumn2,numberColumns_*sizeof(double));
338      delete [] saveRow1;
339      delete [] saveRow2;
340      delete [] saveColumn1;
341      delete [] saveColumn2;
342      objectiveValue_=saveValue;
343    }
344#endif
345#ifdef CLP_DEBUG
346    {
347      int iSequence, number=numberRows_+numberColumns_;
348      for (iSequence=0;iSequence<number;iSequence++) {
349        double lowerValue=lower_[iSequence];
350        double upperValue=upper_[iSequence];
351        double value=solution_[iSequence];
352        if(getStatus(iSequence)!=basic) {
353          assert(lowerValue>-1.0e20);
354          assert(upperValue<1.0e20);
355        }
356        switch(getStatus(iSequence)) {
357         
358        case basic:
359          break;
360        case isFree:
361        case superBasic:
362          break;
363        case atUpperBound:
364          assert (fabs(value-upperValue)<=primalTolerance_) ;
365          break;
366        case atLowerBound:
367          assert (fabs(value-lowerValue)<=primalTolerance_) ;
368          break;
369        }
370      }
371    }
372    if(numberIterations_==debugIteration) {
373      printf("dodgy iteration coming up\n");
374    }
375#endif
376    // choose row to go out
377    // dualRow will go to virtual row pivot choice algorithm
378    dualRow();
379    if (pivotRow_>=0) {
380      // we found a pivot row
381      handler_->message(CLP_SIMPLEX_PIVOTROW,messages_)
382        <<pivotRow_
383        <<CoinMessageEol;
384      // check accuracy of weights
385      dualRowPivot_->checkAccuracy();
386      // get sign for finding row of tableau
387      rowArray_[0]->insert(pivotRow_,directionOut_);
388      factorization_->updateColumnTranspose(rowArray_[1],rowArray_[0]);
389      // put row of tableau in rowArray[0] and columnArray[0]
390      matrix_->transposeTimes(this,-1.0,
391                              rowArray_[0],columnArray_[1],columnArray_[0]);
392      // rowArray has pi equivalent
393      // do ratio test
394      dualColumn(rowArray_[0],columnArray_[0],columnArray_[1],
395                 rowArray_[3]);
396      if (sequenceIn_>=0) {
397        // normal iteration
398        // update the incoming column
399        unpack(rowArray_[1]);
400        factorization_->updateColumn(rowArray_[2],rowArray_[1],true);
401        // and update dual weights (can do in parallel - with extra array)
402        dualRowPivot_->updateWeights(rowArray_[0],rowArray_[2],
403                                     rowArray_[1]);
404        // see if update stable
405        double btranAlpha = -alpha_*directionOut_; // for check
406        alpha_=(*rowArray_[1])[pivotRow_];
407#ifdef CLP_DEBUG
408        if ((handler_->logLevel()&32))
409          printf("btran alpha %g, ftran alpha %g\n",btranAlpha,alpha_);
410#endif
411        double checkValue=1.0e-7;
412        // if can't trust much and long way from optimal then relax
413        if (largestPrimalError_>10.0)
414          checkValue = min(1.0e-4,1.0e-8*largestPrimalError_);
415        if (fabs(btranAlpha)<1.0e-12||fabs(alpha_)<1.0e-12||
416            fabs(btranAlpha-alpha_)>checkValue*(1.0+fabs(alpha_))) {
417          handler_->message(CLP_DUAL_CHECK,messages_)
418            <<btranAlpha
419            <<alpha_
420            <<CoinMessageEol;
421          if (factorization_->pivots()) {
422            dualRowPivot_->unrollWeights();
423            problemStatus_=-2; // factorize now
424            rowArray_[0]->clear();
425            rowArray_[1]->clear();
426            columnArray_[0]->clear();
427            returnCode=-2;
428            break;
429          } else {
430            // take on more relaxed criterion
431            if (fabs(btranAlpha)<1.0e-12||fabs(alpha_)<1.0e-12||
432                fabs(btranAlpha-alpha_)>1.0e-4*(1.0+fabs(alpha_))) {
433              dualRowPivot_->unrollWeights();
434              // need to reject something
435              char x = isColumn(sequenceOut_) ? 'C' :'R';
436              handler_->message(CLP_SIMPLEX_FLAG,messages_)
437                <<x<<sequenceWithin(sequenceOut_)
438                <<CoinMessageEol;
439              setFlagged(sequenceOut_);
440              lastBadIteration_ = numberIterations_; // say be more cautious
441              rowArray_[0]->clear();
442              rowArray_[1]->clear();
443              columnArray_[0]->clear();
444              continue;
445            }
446          }
447        }
448        // update duals BEFORE replaceColumn so can do updateColumn
449        double objectiveChange=0.0;
450        // do duals first as variables may flip bounds
451        // rowArray_[0] and columnArray_[0] may have flips
452        // so use rowArray_[3] for work array from here on
453        int nswapped = 
454          updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[2],theta_,
455                            objectiveChange);
456        // which will change basic solution
457        if (nswapped) {
458#ifdef CLP_DEBUG
459          if ((handler_->logLevel()&16))
460            printf("old dualOut_ %g, v %g, l %g, u %g - new ",
461                   dualOut_,valueOut_,lowerOut_,upperOut_);
462          double oldOut=dualOut_;
463#endif
464          factorization_->updateColumn(rowArray_[3],rowArray_[2],false);
465          dualRowPivot_->updatePrimalSolution(rowArray_[2],
466                                              1.0,objectiveChange);
467         
468          // recompute dualOut_
469          valueOut_ = solution_[sequenceOut_];
470          if (directionOut_<0) {
471            dualOut_ = valueOut_ - upperOut_;
472          } else {
473            dualOut_ = lowerOut_ - valueOut_;
474          }
475#ifdef CLP_DEBUG
476          if ((handler_->logLevel()&16))
477            printf("%g\n",dualOut_);
478          assert(dualOut_<=oldOut);
479#endif
480          if(dualOut_<-10.0e-8&&factorization_->pivots()&&
481             getStatus(sequenceIn_)!=isFree) {
482            // going backwards - factorize
483            dualRowPivot_->unrollWeights();
484            problemStatus_=-2; // factorize now
485            returnCode=-2;
486            break;
487          }
488        }
489        // amount primal will move
490        double movement = -dualOut_*directionOut_/alpha_;
491        // so objective should increase by fabs(dj)*movement
492        // but we already have objective change - so check will be good
493        if (objectiveChange+fabs(movement*dualIn_)<-1.0e-5) {
494#ifdef CLP_DEBUG
495          if (handler_->logLevel()&32)
496            printf("movement %g, swap change %g, rest %g  * %g\n",
497                   objectiveChange+fabs(movement*dualIn_),
498                   objectiveChange,movement,dualIn_);
499#endif
500          if(factorization_->pivots()>5) {
501            // going backwards - factorize
502            dualRowPivot_->unrollWeights();
503            problemStatus_=-2; // factorize now
504            returnCode=-2;
505            break;
506          }
507        }
508        assert(fabs(dualOut_)<1.0e50);
509        // if stable replace in basis
510        int updateStatus = factorization_->replaceColumn(rowArray_[2],
511                                                         pivotRow_,
512                                                         alpha_);
513        // if no pivots, bad update but reasonable alpha - take and invert
514        if (updateStatus==2&&
515                   !factorization_->pivots()&&fabs(alpha_)>1.0e-5)
516          updateStatus=4;
517        if (updateStatus==1||updateStatus==4) {
518          // slight error
519          if (factorization_->pivots()>5||updateStatus==4) {
520            problemStatus_=-2; // factorize now
521            returnCode=-3;
522          }
523        } else if (updateStatus==2) {
524          // major error
525          dualRowPivot_->unrollWeights();
526          // later we may need to unwind more e.g. fake bounds
527          if (factorization_->pivots()) {
528            problemStatus_=-2; // factorize now
529            returnCode=-2;
530            break;
531          } else {
532            // need to reject something
533            char x = isColumn(sequenceOut_) ? 'C' :'R';
534            handler_->message(CLP_SIMPLEX_FLAG,messages_)
535              <<x<<sequenceWithin(sequenceOut_)
536              <<CoinMessageEol;
537            setFlagged(sequenceOut_);
538            lastBadIteration_ = numberIterations_; // say be more cautious
539            rowArray_[0]->clear();
540            rowArray_[1]->clear();
541            columnArray_[0]->clear();
542            // make sure dual feasible
543            // look at all rows and columns
544            CoinIotaN(rowArray_[0]->getIndices(),numberRows_,0);
545            rowArray_[0]->setNumElements(numberRows_);
546            CoinIotaN(columnArray_[0]->getIndices(),numberColumns_,0);
547            columnArray_[0]->setNumElements(numberColumns_);
548            double objectiveChange=0.0;
549            updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[1],
550                              0.0,objectiveChange);
551            continue;
552          }
553        } else if (updateStatus==3) {
554          // out of memory
555          // increase space if not many iterations
556          if (factorization_->pivots()<
557              0.5*factorization_->maximumPivots()&&
558              factorization_->pivots()<200)
559            factorization_->areaFactor(
560                                       factorization_->areaFactor() * 1.1);
561          problemStatus_=-2; // factorize now
562        } 
563        // update primal solution
564        if (theta_<0.0) {
565#ifdef CLP_DEBUG
566          if (handler_->logLevel()&32)
567            printf("negative theta %g\n",theta_);
568#endif
569          theta_=0.0;
570        }
571        // do actual flips
572        flipBounds(rowArray_[0],columnArray_[0],theta_);
573        dualRowPivot_->updatePrimalSolution(rowArray_[1],
574                                            movement,
575                                            objectiveChange);
576#ifdef CLP_DEBUG
577        double oldobj=objectiveValue_;
578#endif
579        // modify dualout
580        dualOut_ /= alpha_;
581        dualOut_ *= -directionOut_;
582        //setStatus(sequenceIn_,basic);
583        dj_[sequenceIn_]=0.0;
584        double oldValue=valueIn_;
585        if (directionIn_==-1) {
586          // as if from upper bound
587          valueIn_ = upperIn_+dualOut_;
588        } else {
589          // as if from lower bound
590          valueIn_ = lowerIn_+dualOut_;
591        }
592        objectiveChange += cost_[sequenceIn_]*(valueIn_-oldValue);
593        // outgoing
594        if (directionOut_>0) {
595          valueOut_ = lowerOut_;
596          //setStatus(sequenceOut_,atLowerBound);
597          dj_[sequenceOut_]=theta_;
598        } else {
599          valueOut_ = upperOut_;
600          //setStatus(sequenceOut_,atUpperBound);
601          dj_[sequenceOut_]=-theta_;
602        }
603        solution_[sequenceOut_]=valueOut_;
604        int whatNext=housekeeping(objectiveChange);
605        // and set bounds correctly
606        originalBound(sequenceIn_); 
607        changeBound(sequenceOut_);
608#ifdef CLP_DEBUG
609        if (objectiveValue_<oldobj-1.0e-5&&(handler_->logLevel()&16))
610          printf("obj backwards %g %g\n",objectiveValue_,oldobj);
611#endif
612        if (whatNext==1) {
613          problemStatus_ =-2; // refactorize
614        } else if (whatNext==2) {
615          // maximum iterations or equivalent
616          problemStatus_= 3;
617          returnCode=3;
618          break;
619        }
620      } else {
621        // no incoming column is valid
622#ifdef CLP_DEBUG
623        if (handler_->logLevel()&32)
624          printf("** no column pivot\n");
625#endif
626        if (factorization_->pivots()<5) {
627          problemStatus_=-4; //say looks infeasible
628          // create ray anyway
629          delete [] ray_;
630          ray_ = new double [ numberRows_];
631          ClpDisjointCopyN(rowArray_[0]->denseVector(),numberRows_,ray_);
632        }
633        rowArray_[0]->clear();
634        columnArray_[0]->clear();
635        returnCode=1;
636        break;
637      }
638    } else {
639      // no pivot row
640#ifdef CLP_DEBUG
641      if (handler_->logLevel()&32)
642        printf("** no row pivot\n");
643#endif
644      if (!factorization_->pivots()) {
645        // may have crept through - so may be optimal
646        // check any flagged variables
647        int iRow;
648        for (iRow=0;iRow<numberRows_;iRow++) {
649          int iPivot=pivotVariable_[iRow];
650          if (flagged(iPivot))
651            break;
652        }
653        if (numberFake_||numberDualInfeasibilities_) {
654          // may be dual infeasible
655          problemStatus_=-5;
656        } else {
657          if (iRow<numberRows_) {
658#ifdef CLP_DEBUG
659            std::cerr<<"Flagged variables at end - infeasible?"<<std::endl;
660            printf("Probably infeasible - pivot was %g\n",alpha_);
661#endif
662            if (fabs(alpha_)<1.0e-4) {
663              problemStatus_=1;
664            } else {
665#ifdef CLP_DEBUG
666              abort();
667#endif
668            }
669          } else {
670            problemStatus_=0;
671          }
672        }
673      } else {
674        problemStatus_=-3;
675      }
676      returnCode=0;
677      break;
678    }
679  }
680  return returnCode;
681}
682/* The duals are updated by the given arrays.
683   Returns number of infeasibilities.
684   rowArray and columnarray will have flipped
685   The output vector has movement (row length array) */
686int
687ClpSimplexDual::updateDualsInDual(CoinIndexedVector * rowArray,
688                                  CoinIndexedVector * columnArray,
689                                  CoinIndexedVector * outputArray,
690                                  double theta,
691                                  double & objectiveChange)
692{
693 
694  outputArray->clear();
695 
696  double * work;
697  int number;
698  int * which;
699 
700  int numberInfeasibilities=0;
701  int numberRowInfeasibilities=0;
702 
703  // see whether we will be doing full recompute
704  bool fullRecompute= (rowArray->getNumElements()==numberRows_&&
705                       columnArray->getNumElements()==numberColumns_);
706  int numberAtFake=0;
707 
708  // use a tighter tolerance except for all being okay
709  double tolerance = dualTolerance_;
710 
711  double changeObj=0.0;
712 
713  int iSection;
714 
715  for (iSection=0;iSection<2;iSection++) {
716    int i;
717    double * solution = solutionRegion(iSection);
718    double * reducedCost = djRegion(iSection);
719    double * lower = lowerRegion(iSection);
720    double * upper = upperRegion(iSection);
721    double * cost = costRegion(iSection);
722    int addSequence;
723    if (!iSection) {
724      addSequence = numberColumns_;
725      work = rowArray->denseVector();
726      number = rowArray->getNumElements();
727      which = rowArray->getIndices();
728    } else {
729      // set number of infeasibilities in row array
730      addSequence=0;
731      numberRowInfeasibilities=numberInfeasibilities;
732      rowArray->setNumElements(numberInfeasibilities);
733      numberInfeasibilities=0;
734      work = columnArray->denseVector();
735      number = columnArray->getNumElements();
736      which = columnArray->getIndices();
737    }
738   
739    for (i=0;i<number;i++) {
740      int iSequence = which[i];
741      double alphaI = work[iSequence];
742      double value = reducedCost[iSequence]-theta*alphaI;
743      work[iSequence]=0.0;
744      reducedCost[iSequence]=value;
745     
746      if (!fixed(iSequence+addSequence)) {
747        double movement=0.0;
748        FakeBound bound = getFakeBound(iSequence+addSequence);
749        Status status = getStatus(iSequence+addSequence);
750
751        switch(status) {
752         
753        case basic:
754        case superBasic:
755          break;
756        case isFree:
757          if (fabs(value)>tolerance) { 
758#ifdef CLP_DEBUG
759            if (handler_->logLevel()&32)
760              printf("%d %d, free has dj of %g, alpha %g\n",
761                     iSection,iSequence,value,alphaI);
762#endif
763          }
764          break;
765        case atUpperBound:
766          if (value>tolerance) {
767            // to lower bound (if swap)
768            // put back alpha
769            which[numberInfeasibilities++]=iSequence;
770            work[iSequence]=alphaI;
771            movement = lower[iSequence]-upper[iSequence];
772#ifdef CLP_DEBUG
773            if ((handler_->logLevel()&32))
774              printf("%d %d, new dj %g, alpha %g, movement %g\n",
775                   iSection,iSequence,value,alphaI,movement);
776#endif
777            changeObj += movement*cost[iSequence];
778            if (bound==ClpSimplexDual::bothFake||
779                bound==ClpSimplexDual::lowerFake) 
780              numberAtFake++;
781          } else if (fullRecompute) {
782            // at correct bound
783            if (bound==ClpSimplexDual::bothFake||
784                bound==ClpSimplexDual::upperFake) {
785              // but flip if dj would allow
786              if (bound==ClpSimplexDual::upperFake&&
787                  value>=-tolerance) {
788                movement = lower[iSequence]-upper[iSequence];
789                setStatus(iSequence+addSequence,atLowerBound);
790                solution[iSequence] = lower[iSequence];
791                changeObj += movement*cost[iSequence];
792              } else {
793                numberAtFake++;
794              }
795            }
796          }
797          break;
798        case atLowerBound:
799          if (value<-tolerance) {
800            // to upper bound
801            // put back alpha
802            which[numberInfeasibilities++]=iSequence;
803            work[iSequence]=alphaI;
804            movement = upper[iSequence] - lower[iSequence];
805#ifdef CLP_DEBUG
806            if ((handler_->logLevel()&32))
807              printf("%d %d, new dj %g, alpha %g, movement %g\n",
808                   iSection,iSequence,value,alphaI,movement);
809#endif
810            changeObj += movement*cost[iSequence];
811            if (bound==ClpSimplexDual::bothFake||
812                bound==ClpSimplexDual::upperFake) 
813              numberAtFake++;
814          } else if (fullRecompute) {
815            // at correct bound
816            if (bound==ClpSimplexDual::bothFake||
817                bound==ClpSimplexDual::lowerFake) {
818              // but flip if dj would allow
819              if (bound==ClpSimplexDual::lowerFake&&
820                  value<=tolerance) {
821                movement = upper[iSequence] - lower[iSequence];
822                setStatus(iSequence+addSequence,atUpperBound);
823                solution[iSequence] = upper[iSequence];
824                changeObj += movement*cost[iSequence];
825              } else {
826                numberAtFake++;
827              }
828            }
829          }
830          break;
831        }
832        if (!fullRecompute) {
833          if (movement) {
834            if (!iSection) {
835              // row (sign ?)
836              outputArray->quickAdd(iSequence,-movement);
837            } else {
838              matrix_->add(this,outputArray,iSequence,movement);
839            }
840          }
841        }
842      }
843    }
844  }
845#ifdef CLP_DEBUG
846  if (fullRecompute&&numberAtFake&&(handler_->logLevel()&16)!=0) 
847    printf("%d fake after full update\n",numberAtFake);
848#endif
849  // set number of infeasibilities
850  columnArray->setNumElements(numberInfeasibilities);
851  numberInfeasibilities += numberRowInfeasibilities;
852  if (fullRecompute) {
853    // do actual flips
854    flipBounds(rowArray,columnArray,theta);
855    numberFake_ = numberAtFake;
856  }
857  objectiveChange += changeObj;
858  return numberInfeasibilities;
859}
860/*
861   Chooses dual pivot row
862   Would be faster with separate region to scan
863   and will have this (with square of infeasibility) when steepest
864   For easy problems we can just choose one of the first rows we look at
865*/
866void
867ClpSimplexDual::dualRow()
868{
869  // get pivot row using whichever method it is
870  // first see if any free variables and put them in basis
871  int nextFree = nextSuperBasic();
872  int chosenRow=-1;
873  //nextFree=-1; //off
874  if (nextFree>=0) {
875    // unpack vector and find a good pivot
876    unpack(rowArray_[1],nextFree);
877    factorization_->updateColumn(rowArray_[2],rowArray_[1],false);
878
879    double * work=rowArray_[1]->denseVector();
880    int number=rowArray_[1]->getNumElements();
881    int * which=rowArray_[1]->getIndices();
882    double bestFeasibleAlpha=0.0;
883    int bestFeasibleRow=-1;
884    double bestInfeasibleAlpha=0.0;
885    int bestInfeasibleRow=-1;
886    int i;
887
888    for (i=0;i<number;i++) {
889      int iRow = which[i];
890      double alpha = fabs(work[iRow]);
891      if (alpha>1.0e-3) {
892        int iSequence=pivotVariable_[iRow];
893        double value = solution_[iSequence];
894        double lower = lower_[iSequence];
895        double upper = upper_[iSequence];
896        double infeasibility=0.0;
897        if (value>upper)
898          infeasibility = value-upper;
899        else if (value<lower)
900          infeasibility = lower-value;
901        if (infeasibility*alpha>bestInfeasibleAlpha&&alpha>1.0e-1) {
902          if (!flagged(iSequence)) {
903            bestInfeasibleAlpha = infeasibility*alpha;
904            bestInfeasibleRow=iRow;
905          }
906        }
907        if (alpha>bestFeasibleAlpha&&(lower>-1.0e20||upper<1.0e20)) {
908          bestFeasibleAlpha = alpha;
909          bestFeasibleRow=iRow;
910        }
911      }
912    }
913    if (bestInfeasibleRow>=0) 
914      chosenRow = bestInfeasibleRow;
915    else if (bestFeasibleAlpha>1.0e-2)
916      chosenRow = bestFeasibleRow;
917    if (chosenRow>=0)
918      pivotRow_=chosenRow;
919    rowArray_[1]->clear();
920  } 
921  if (chosenRow<0) 
922    pivotRow_=dualRowPivot_->pivotRow();
923
924  if (pivotRow_>=0) {
925    sequenceOut_ = pivotVariable_[pivotRow_];
926    valueOut_ = solution_[sequenceOut_];
927    lowerOut_ = lower_[sequenceOut_];
928    upperOut_ = upper_[sequenceOut_];
929    // if we have problems we could try other way and hope we get a
930    // zero pivot?
931    if (valueOut_>upperOut_) {
932      directionOut_ = -1;
933      dualOut_ = valueOut_ - upperOut_;
934    } else if (valueOut_<lowerOut_) {
935      directionOut_ = 1;
936      dualOut_ = lowerOut_ - valueOut_;
937    } else {
938      // odd (could be free) - it's feasible - go to nearest
939      if (valueOut_-lowerOut_<upperOut_-valueOut_) {
940        directionOut_ = 1;
941        dualOut_ = lowerOut_ - valueOut_;
942      } else {
943        directionOut_ = -1;
944        dualOut_ = valueOut_ - upperOut_;
945      }
946    }
947#ifdef CLP_DEBUG
948    assert(dualOut_>=0.0);
949#endif
950  }
951  return ;
952}
953// Checks if any fake bounds active - if so returns number and modifies
954// dualBound_ and everything.
955// Free variables will be left as free
956// Returns number of bounds changed if >=0
957// Returns -1 if not initialize and no effect
958// Fills in changeVector which can be used to see if unbounded
959// and cost of change vector
960int
961ClpSimplexDual::changeBounds(bool initialize,
962                                 CoinIndexedVector * outputArray,
963                                 double & changeCost)
964{
965  if (!initialize) {
966    int numberInfeasibilities;
967    double newBound;
968    newBound = 5.0*dualBound_;
969    numberInfeasibilities=0;
970    changeCost=0.0;
971    // put back original bounds and then check
972    createRim(3);
973    int iSequence;
974    // bounds will get bigger - just look at ones at bounds
975    for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
976      double lowerValue=lower_[iSequence];
977      double upperValue=upper_[iSequence];
978      double value=solution_[iSequence];
979      setFakeBound(iSequence,ClpSimplexDual::noFake);
980      switch(getStatus(iSequence)) {
981
982      case basic:
983        break;
984      case isFree:
985      case superBasic:
986        break;
987      case atUpperBound:
988        if (fabs(value-upperValue)>primalTolerance_) 
989          numberInfeasibilities++;
990        break;
991      case atLowerBound:
992        if (fabs(value-lowerValue)>primalTolerance_) 
993          numberInfeasibilities++;
994        break;
995      }
996    }
997    if (numberInfeasibilities) {
998      int iSequence;
999      for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
1000        double lowerValue=lower_[iSequence];
1001        double upperValue=upper_[iSequence];
1002        double newLowerValue;
1003        double newUpperValue;
1004        Status status = getStatus(iSequence);
1005        if (status==atUpperBound||
1006            status==atLowerBound) {
1007          double value = solution_[iSequence];
1008          if (value-lowerValue<=upperValue-value) {
1009            newLowerValue = max(lowerValue,value-0.666667*newBound);
1010            newUpperValue = min(upperValue,newLowerValue+newBound);
1011          } else {
1012            newUpperValue = min(upperValue,value+0.666667*newBound);
1013            newLowerValue = max(lowerValue,newUpperValue-newBound);
1014          }
1015          lower_[iSequence]=newLowerValue;
1016          upper_[iSequence]=newUpperValue;
1017          if (newLowerValue > lowerValue) {
1018            if (newUpperValue < upperValue) {
1019              setFakeBound(iSequence,ClpSimplexDual::bothFake);
1020            } else {
1021              setFakeBound(iSequence,ClpSimplexDual::lowerFake);
1022            }
1023          } else {
1024            if (newUpperValue < upperValue) 
1025              setFakeBound(iSequence,ClpSimplexDual::upperFake);
1026          }
1027          if (status==atUpperBound)
1028            solution_[iSequence] = newUpperValue;
1029          else
1030            solution_[iSequence] = newLowerValue;
1031          double movement = solution_[iSequence] - value;
1032          if (movement&&outputArray) {
1033            if (iSequence>=numberColumns_) {
1034              outputArray->quickAdd(iSequence,-movement);
1035              changeCost += movement*cost_[iSequence];
1036            } else {
1037              matrix_->add(this,outputArray,iSequence,movement);
1038              changeCost += movement*cost_[iSequence];
1039            }
1040          }
1041        }
1042      }
1043      dualBound_ = newBound;
1044    } else {
1045      numberInfeasibilities=-1;
1046    }
1047    return numberInfeasibilities;
1048  } else {
1049    int iSequence;
1050     
1051    for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
1052      Status status = getStatus(iSequence);
1053      if (status==atUpperBound||
1054          status==atLowerBound) {
1055        double lowerValue=lower_[iSequence];
1056        double upperValue=upper_[iSequence];
1057        double value = solution_[iSequence];
1058        if (lowerValue>-largeValue_||upperValue<largeValue_) {
1059          if (lowerValue-value>-0.5*dualBound_||
1060              upperValue-value<0.5*dualBound_) {
1061            if (fabs(lowerValue-value)<=fabs(upperValue-value)) {
1062              if (upperValue > lowerValue + dualBound_) {
1063                upper_[iSequence]=lowerValue+dualBound_;
1064                setFakeBound(iSequence,ClpSimplexDual::upperFake);
1065              }
1066            } else {
1067              if (lowerValue < upperValue - dualBound_) {
1068                lower_[iSequence]=upperValue-dualBound_;
1069                setFakeBound(iSequence,ClpSimplexDual::lowerFake);
1070              }
1071            }
1072          } else {
1073            lower_[iSequence]=-0.5*dualBound_;
1074            upper_[iSequence]= 0.5*dualBound_;
1075            setFakeBound(iSequence,ClpSimplexDual::bothFake);
1076          }
1077        } else {
1078          // set non basic free variables to fake bounds
1079          // I don't think we should ever get here
1080          assert(!("should not be here"));
1081          lower_[iSequence]=-0.5*dualBound_;
1082          upper_[iSequence]= 0.5*dualBound_;
1083          setFakeBound(iSequence,ClpSimplexDual::bothFake);
1084          setStatus(iSequence,atUpperBound);
1085          solution_[iSequence]=0.5*dualBound_;
1086        }
1087      }
1088    }
1089
1090    return 1;
1091  }
1092}
1093/*
1094   Row array has row part of pivot row (as duals so sign may be switched)
1095   Column array has column part.
1096   This chooses pivot column.
1097   Spare array will be needed when we start getting clever.
1098   We will check for basic so spare array will never overflow.
1099   If necessary will modify costs
1100*/
1101void
1102ClpSimplexDual::dualColumn(CoinIndexedVector * rowArray,
1103                       CoinIndexedVector * columnArray,
1104                       CoinIndexedVector * spareArray,
1105                       CoinIndexedVector * spareArray2)
1106{
1107  double * work;
1108  int number;
1109  int * which;
1110  double * reducedCost;
1111
1112  int iSection;
1113
1114  sequenceIn_=-1;
1115  int numberPossiblySwapped=0;
1116  int numberRemaining=0;
1117 
1118  double totalThru=0.0; // for when variables flip
1119  double acceptablePivot=1.0e-7;
1120  if (factorization_->pivots()>10)
1121    acceptablePivot=1.0e-5; // if we have iterated be more strict
1122  else if (factorization_->pivots()>5)
1123    acceptablePivot=1.0e-6; // if we have iterated be slightly more strict
1124  double bestEverPivot=acceptablePivot;
1125  int lastSequence = -1;
1126  double lastPivot=0.0;
1127  double upperTheta;
1128  double newTolerance = dualTolerance_;
1129  // will we need to increase tolerance
1130  bool thisIncrease=false;
1131  // If we think we need to modify costs (not if something from broad sweep)
1132  bool modifyCosts=false;
1133  // Increase in objective due to swapping bounds (may be negative)
1134  double increaseInObjective=0.0;
1135
1136  // use spareArrays to put ones looked at in
1137  // we are going to flip flop between
1138  int iFlip = 0;
1139  // Possible list of pivots
1140  int interesting[2];
1141  // where possible swapped ones are
1142  int swapped[2];
1143  // for zeroing out arrays after
1144  int marker[2][2];
1145  // pivot elements
1146  double * array[2], * spare, * spare2;
1147  // indices
1148  int * indices[2], * index, * index2;
1149  spareArray->clear();
1150  spareArray2->clear();
1151  array[0] = spareArray->denseVector();
1152  indices[0] = spareArray->getIndices();
1153  spare = array[0];
1154  index = indices[0];
1155  array[1] = spareArray2->denseVector();
1156  indices[1] = spareArray2->getIndices();
1157  int i;
1158  double * lower;
1159  double * upper;
1160
1161  // initialize lists
1162  for (i=0;i<2;i++) {
1163    interesting[i]=0;
1164    swapped[i]=numberColumns_;
1165    marker[i][0]=0;
1166    marker[i][1]=numberColumns_;
1167  }
1168
1169  /*
1170    First we get a list of possible pivots.  We can also see if the
1171    problem looks infeasible or whether we want to pivot in free variable.
1172    This may make objective go backwards but can only happen a finite
1173    number of times and I do want free variables basic.
1174
1175    Then we flip back and forth.  At the start of each iteration
1176    interesting[iFlip] should have possible candidates and swapped[iFlip]
1177    will have pivots if we decide to take a previous pivot.
1178    At end of each iteration interesting[1-iFlip] should have
1179    candidates if we go through this theta and swapped[1-iFlip]
1180    pivots if we don't go through.
1181
1182    At first we increase theta and see what happens.  We start
1183    theta at a reasonable guess.  If in right area then we do bit by bit.
1184
1185   */
1186
1187  // do first pass to get possibles
1188  // We can also see if infeasible or pivoting on free
1189  double tentativeTheta = 1.0e22;
1190  upperTheta = 1.0e31;
1191  double freePivot = acceptablePivot;
1192  for (iSection=0;iSection<2;iSection++) {
1193
1194    int addSequence;
1195
1196    if (!iSection) {
1197      lower = rowLowerWork_;
1198      upper = rowUpperWork_;
1199      work = rowArray->denseVector();
1200      number = rowArray->getNumElements();
1201      which = rowArray->getIndices();
1202      reducedCost = rowReducedCost_;
1203      addSequence = numberColumns_;
1204    } else {
1205      lower = columnLowerWork_;
1206      upper = columnUpperWork_;
1207      work = columnArray->denseVector();
1208      number = columnArray->getNumElements();
1209      which = columnArray->getIndices();
1210      reducedCost = reducedCostWork_;
1211      addSequence = 0;
1212    }
1213   
1214    for (i=0;i<number;i++) {
1215      int iSequence = which[i];
1216      double alpha = work[iSequence];
1217      if (fixed(iSequence+addSequence)||!alpha) 
1218        continue; // skip fixed ones or (zeroed out)
1219      double oldValue = reducedCost[iSequence];
1220      double value = oldValue-tentativeTheta*alpha;
1221      int keep = 0;
1222
1223      switch(getStatus(iSequence+addSequence)) {
1224         
1225      case basic:
1226        break;
1227      case isFree:
1228      case superBasic:
1229        if (oldValue>dualTolerance_) {
1230          keep = 2;
1231        } else if (oldValue<-dualTolerance_) {
1232          keep = 2;
1233        } else {
1234          if (fabs(alpha)>10.0*acceptablePivot) 
1235            keep = 2;
1236
1237
1238        }
1239        break;
1240      case atUpperBound:
1241        assert (oldValue<=dualTolerance_*1.0001);
1242        if (value>newTolerance) {
1243          keep = 1;
1244          value = oldValue-upperTheta*alpha;
1245          if (value>newTolerance && -alpha>=acceptablePivot) 
1246            upperTheta = (oldValue-newTolerance)/alpha;
1247        }
1248        break;
1249      case atLowerBound:
1250        assert (oldValue>=-dualTolerance_*1.0001);
1251        if (value<-newTolerance) {
1252          keep = 1;
1253          value = oldValue-upperTheta*alpha;
1254          if (value<-newTolerance && alpha>=acceptablePivot) 
1255            upperTheta = (oldValue+newTolerance)/alpha;
1256        }
1257        break;
1258      }
1259      if (keep) {
1260        if (keep==2) {
1261          // free - choose largest
1262          if (fabs(alpha)>freePivot) {
1263            freePivot=fabs(alpha);
1264            sequenceIn_ = iSequence + addSequence;
1265            theta_=oldValue/alpha;
1266          }
1267        } else {
1268          // add to list
1269          spare[numberRemaining]=alpha;
1270          index[numberRemaining++]=iSequence+addSequence;
1271        }
1272      }
1273    }
1274  }
1275  interesting[0]=numberRemaining;
1276  marker[0][0] = numberRemaining;
1277
1278  if (!numberRemaining&&sequenceIn_<0)
1279    return; // Looks infeasible
1280
1281  if (sequenceIn_>=0) {
1282    // free variable - always choose
1283  } else {
1284
1285    theta_=1.0e50;
1286    // now flip flop between spare arrays until reasonable theta
1287    tentativeTheta = max(10.0*upperTheta,1.0e-7);
1288   
1289    // loops increasing tentative theta until can't go through
1290   
1291    while (tentativeTheta < 1.0e22) {
1292      double thruThis = 0.0;
1293     
1294      double bestPivot=acceptablePivot;
1295      int bestSequence=-1;
1296     
1297      numberPossiblySwapped = numberColumns_;
1298      numberRemaining = 0;
1299
1300      upperTheta = 1.0e50;
1301
1302      spare = array[iFlip];
1303      index = indices[iFlip];
1304      spare2 = array[1-iFlip];
1305      index2 = indices[1-iFlip];
1306
1307      // try 3 different ways
1308      // 1 bias increase by ones with slightly wrong djs
1309      // 2 bias by all
1310      // 3 bias by all - tolerance (doesn't seem very good)
1311#define TRYBIAS 1
1312
1313
1314      double increaseInThis=0.0; //objective increase in this loop
1315     
1316      for (i=0;i<interesting[iFlip];i++) {
1317        int iSequence = index[i];
1318        double alpha = spare[i];
1319        double oldValue = dj_[iSequence];
1320        double value = oldValue-tentativeTheta*alpha;
1321
1322        if (alpha < 0.0) {
1323          //at upper bound
1324          if (value>newTolerance) {
1325            double range = upper_[iSequence] - lower_[iSequence];
1326            thruThis -= range*alpha;
1327#if TRYBIAS==1
1328            if (oldValue>0.0)
1329              increaseInThis -= oldValue*range;
1330#elif TRYBIAS==2
1331            increaseInThis -= oldValue*range;
1332#else
1333            increaseInThis -= (oldValue+dualTolerance_)*range;
1334#endif
1335            // goes on swapped list (also means candidates if too many)
1336            spare2[--numberPossiblySwapped]=alpha;
1337            index2[numberPossiblySwapped]=iSequence;
1338            if (fabs(alpha)>bestPivot) {
1339              bestPivot=fabs(alpha);
1340              bestSequence=numberPossiblySwapped;
1341            }
1342          } else {
1343            value = oldValue-upperTheta*alpha;
1344            if (value>newTolerance && -alpha>=acceptablePivot) 
1345              upperTheta = (oldValue-newTolerance)/alpha;
1346            spare2[numberRemaining]=alpha;
1347            index2[numberRemaining++]=iSequence;
1348          }
1349        } else {
1350          // at lower bound
1351          if (value<-newTolerance) {
1352            double range = upper_[iSequence] - lower_[iSequence];
1353            thruThis += range*alpha;
1354            //?? is this correct - and should we look at good ones
1355#if TRYBIAS==1
1356            if (oldValue<0.0)
1357              increaseInThis += oldValue*range;
1358#elif TRYBIAS==2
1359            increaseInThis += oldValue*range;
1360#else
1361            increaseInThis += (oldValue-dualTolerance_)*range;
1362#endif
1363            // goes on swapped list (also means candidates if too many)
1364            spare2[--numberPossiblySwapped]=alpha;
1365            index2[numberPossiblySwapped]=iSequence;
1366            if (fabs(alpha)>bestPivot) {
1367              bestPivot=fabs(alpha);
1368              bestSequence=numberPossiblySwapped;
1369            }
1370          } else {
1371            value = oldValue-upperTheta*alpha;
1372            if (value<-newTolerance && alpha>=acceptablePivot) 
1373              upperTheta = (oldValue+newTolerance)/alpha;
1374            spare2[numberRemaining]=alpha;
1375            index2[numberRemaining++]=iSequence;
1376          }
1377        }
1378      }
1379      swapped[1-iFlip]=numberPossiblySwapped;
1380      interesting[1-iFlip]=numberRemaining;
1381      marker[1-iFlip][0]= max(marker[1-iFlip][0],numberRemaining);
1382      marker[1-iFlip][1]= min(marker[1-iFlip][1],numberPossiblySwapped);
1383     
1384      if (totalThru+thruThis>=fabs(dualOut_)||
1385          increaseInObjective+increaseInThis<0.0) {
1386        // We should be pivoting in this batch
1387        // so compress down to this lot
1388        numberRemaining=0;
1389        for (i=numberColumns_-1;i>=swapped[1-iFlip];i--) {
1390          spare[numberRemaining]=spare2[i];
1391          index[numberRemaining++]=index2[i];
1392        }
1393        interesting[iFlip]=numberRemaining;
1394        int iTry;
1395#define MAXTRY 100
1396        // first get ratio with tolerance
1397        for (iTry=0;iTry<MAXTRY;iTry++) {
1398         
1399          upperTheta=1.0e50;
1400          numberPossiblySwapped = numberColumns_;
1401          numberRemaining = 0;
1402
1403          increaseInThis=0.0; //objective increase in this loop
1404
1405          thruThis=0.0;
1406
1407          spare = array[iFlip];
1408          index = indices[iFlip];
1409          spare2 = array[1-iFlip];
1410          index2 = indices[1-iFlip];
1411     
1412          for (i=0;i<interesting[iFlip];i++) {
1413            int iSequence=index[i];
1414            double alpha=spare[i];
1415            double oldValue = dj_[iSequence];
1416            double value = oldValue-upperTheta*alpha;
1417           
1418            if (alpha < 0.0) {
1419              //at upper bound
1420              if (value>newTolerance) {
1421                if (-alpha>=acceptablePivot) {
1422                  upperTheta = (oldValue-newTolerance)/alpha;
1423                  // recompute value and make sure works
1424                  value = oldValue-upperTheta*alpha;
1425                  if (value<0)
1426                    upperTheta *= 1.0 +1.0e-11; // must be large
1427                }
1428              }
1429            } else {
1430              // at lower bound
1431              if (value<-newTolerance) {
1432                if (alpha>=acceptablePivot) {
1433                  upperTheta = (oldValue+newTolerance)/alpha;
1434                  // recompute value and make sure works
1435                  value = oldValue-upperTheta*alpha;
1436                  if (value>0)
1437                    upperTheta *= 1.0 +1.0e-11; // must be large
1438                }
1439              }
1440            }
1441          }
1442          bestPivot=acceptablePivot;
1443          sequenceIn_=-1;
1444          // now choose largest and sum all ones which will go through
1445#define MINIMUMTHETA 1.0e-12
1446          for (i=0;i<interesting[iFlip];i++) {
1447            int iSequence=index[i];
1448            double alpha=spare[i];
1449            double value = dj_[iSequence]-upperTheta*alpha;
1450            double badDj=0.0;
1451           
1452            bool addToSwapped=false;
1453           
1454            if (alpha < 0.0) {
1455              //at upper bound
1456              if (value>=0.0) { 
1457                addToSwapped=true;
1458#if TRYBIAS==1
1459                badDj = -max(dj_[iSequence],0.0);
1460#elif TRYBIAS==2
1461                badDj = -dj_[iSequence];
1462#else
1463                badDj = -dj_[iSequence]-dualTolerance_;
1464#endif
1465              }
1466            } else {
1467              // at lower bound
1468              if (value<=0.0) {
1469                addToSwapped=true;
1470#if TRYBIAS==1
1471                badDj = min(dj_[iSequence],0.0);
1472#elif TRYBIAS==2
1473                badDj = dj_[iSequence];
1474#else
1475                badDj = dj_[iSequence]-dualTolerance_;
1476#endif
1477              }
1478            }
1479            if (!addToSwapped) {
1480              // add to list of remaining
1481              spare2[numberRemaining]=alpha;
1482              index2[numberRemaining++]=iSequence;
1483            } else {
1484              // add to list of swapped
1485              spare2[--numberPossiblySwapped]=alpha;
1486              index2[numberPossiblySwapped]=iSequence;
1487              // select if largest pivot
1488              if (fabs(alpha)>bestPivot) {
1489                sequenceIn_ = numberPossiblySwapped;
1490                bestPivot =  fabs(alpha);
1491                theta_ = dj_[iSequence]/alpha;
1492              }
1493              double range = upper[iSequence] - lower[iSequence];
1494              thruThis += range*fabs(alpha);
1495              increaseInThis += badDj*range;
1496            }
1497          }
1498          swapped[1-iFlip]=numberPossiblySwapped;
1499          interesting[1-iFlip]=numberRemaining;
1500          marker[1-iFlip][0]= max(marker[1-iFlip][0],numberRemaining);
1501          marker[1-iFlip][1]= min(marker[1-iFlip][1],numberPossiblySwapped);
1502          // If we stop now this will be increase in objective (I think)
1503          double increase = (fabs(dualOut_)-totalThru)*theta_;
1504          increase += increaseInObjective;
1505          if (theta_<0.0)
1506            thruThis += fabs(dualOut_); // force using this one
1507          if (increaseInObjective<0.0&&increase<0.0&&lastSequence>=0) {
1508            // back
1509            // We may need to be more careful - we could do by
1510            // switch so we always do fine grained?
1511            bestPivot=0.0;
1512          } else {
1513            // add in
1514            totalThru += thruThis;
1515            increaseInObjective += increaseInThis;
1516          }
1517          if (bestPivot<0.1*bestEverPivot&&
1518              bestEverPivot>1.0e-6&&bestPivot<1.0e-3) {
1519            // back to previous one
1520            sequenceIn_=lastSequence;
1521            // swap regions
1522            iFlip = 1-iFlip;
1523            break;
1524          } else if (sequenceIn_==-1&&upperTheta>largeValue_) {
1525            if (lastPivot>acceptablePivot) {
1526              // back to previous one
1527              sequenceIn_=lastSequence;
1528              // swap regions
1529              iFlip = 1-iFlip;
1530            } else {
1531              // can only get here if all pivots too small
1532            }
1533            break;
1534          } else if (totalThru>=fabs(dualOut_)) {
1535            modifyCosts=true; // fine grain - we can modify costs
1536            break; // no point trying another loop
1537          } else {
1538            lastSequence=sequenceIn_;
1539            if (bestPivot>bestEverPivot)
1540              bestEverPivot=bestPivot;
1541            iFlip = 1 -iFlip;
1542            modifyCosts=true; // fine grain - we can modify costs
1543        }
1544        }
1545        if (iTry==MAXTRY)
1546          iFlip = 1-iFlip; // flip back
1547        break;
1548      } else {
1549        // skip this lot
1550        if (bestPivot>1.0e-3||bestPivot>bestEverPivot) {
1551          bestEverPivot=bestPivot;
1552          lastSequence=bestSequence;
1553        } else {
1554          // keep old swapped
1555          memcpy(array[1-iFlip]+swapped[iFlip],
1556                 array[iFlip]+swapped[iFlip],
1557                 (numberColumns_-swapped[iFlip])*sizeof(double));
1558          memcpy(indices[1-iFlip]+swapped[iFlip],
1559                 indices[iFlip]+swapped[iFlip],
1560                 (numberColumns_-swapped[iFlip])*sizeof(int));
1561          marker[1-iFlip][1] = min(marker[1-iFlip][1],swapped[iFlip]);
1562          swapped[1-iFlip]=swapped[iFlip];
1563        }
1564        increaseInObjective += increaseInThis;
1565        iFlip = 1 - iFlip; // swap regions
1566        tentativeTheta = 2.0*upperTheta;
1567        totalThru += thruThis;
1568      }
1569    }
1570   
1571    // can get here without sequenceIn_ set but with lastSequence
1572    if (sequenceIn_<0&&lastSequence>=0) {
1573      // back to previous one
1574      sequenceIn_=lastSequence;
1575      // swap regions
1576      iFlip = 1-iFlip;
1577    }
1578   
1579    if (sequenceIn_>=0) {
1580      // at this stage sequenceIn_ is just pointer into index array
1581      // flip just so we can use iFlip
1582      iFlip = 1 -iFlip;
1583      spare = array[iFlip];
1584      index = indices[iFlip];
1585      double oldValue;
1586      double alpha = spare[sequenceIn_];
1587      sequenceIn_ = indices[iFlip][sequenceIn_];
1588      oldValue = dj_[sequenceIn_];
1589      theta_ = oldValue/alpha;
1590      if (theta_<MINIMUMTHETA) {
1591        // can't pivot to zero
1592        if (oldValue-MINIMUMTHETA*alpha>=-dualTolerance_) {
1593          theta_=MINIMUMTHETA;
1594        } else if (oldValue-MINIMUMTHETA*alpha>=-newTolerance) {
1595          theta_=MINIMUMTHETA;
1596          thisIncrease=true;
1597        } else {
1598          theta_=(oldValue+newTolerance)/alpha;
1599          assert(theta_>=0.0);
1600          thisIncrease=true;
1601        } 
1602      }
1603      // may need to adjust costs so all dual feasible AND pivoted is exactly 0
1604      if (modifyCosts) {
1605        int i;
1606        double * workRow = rowArray->denseVector();
1607        double * workColumn = columnArray->denseVector();
1608        for (i=numberColumns_-1;i>=swapped[iFlip];i--) {
1609          int iSequence=index[i];
1610          double alpha;
1611          if (iSequence>=numberColumns_)
1612            alpha=workRow[iSequence-numberColumns_];
1613          else
1614            alpha=workColumn[iSequence];
1615          double value = dj_[iSequence]-theta_*alpha;
1616           
1617          // can't be free here
1618         
1619          if (alpha < 0.0) {
1620            //at upper bound
1621            if (value>dualTolerance_) {
1622              thisIncrease=true;
1623#define MODIFYCOST 2
1624#if MODIFYCOST
1625              // modify cost to hit new tolerance
1626              double modification = alpha*theta_-dj_[iSequence]
1627                +newTolerance;
1628              //modification = min(modification,dualTolerance_);
1629              //assert (fabs(modification)<1.0e-7);
1630              dj_[iSequence] += modification;
1631              cost_[iSequence] +=  modification;
1632#endif
1633            }
1634          } else {
1635            // at lower bound
1636            if (-value>dualTolerance_) {
1637              thisIncrease=true;
1638#if MODIFYCOST
1639              // modify cost to hit new tolerance
1640              double modification = alpha*theta_-dj_[iSequence]
1641                -newTolerance;
1642              //modification = max(modification,-dualTolerance_);
1643              //assert (fabs(modification)<1.0e-7);
1644              dj_[iSequence] += modification;
1645              cost_[iSequence] +=  modification;
1646#endif
1647            }
1648          }
1649        }
1650      }
1651    }
1652  }
1653
1654  if (sequenceIn_>=0) {
1655    if (sequenceIn_>=numberColumns_) {
1656      //slack
1657      alpha_ = rowArray->denseVector()[sequenceIn_-numberColumns_];
1658    } else {
1659      // column
1660      alpha_ = columnArray->denseVector()[sequenceIn_];
1661    }
1662    lowerIn_ = lower_[sequenceIn_];
1663    upperIn_ = upper_[sequenceIn_];
1664    valueIn_ = solution_[sequenceIn_];
1665    dualIn_ = dj_[sequenceIn_];
1666
1667    if (numberTimesOptimal_) {
1668      // can we adjust cost back closer to original
1669      //*** add coding
1670    }
1671#if MODIFYCOST>1   
1672    // modify cost to hit zero exactly
1673    // so (dualIn_+modification)==theta_*alpha_
1674    double modification = theta_*alpha_-dualIn_;
1675    dualIn_ += modification;
1676    dj_[sequenceIn_]=dualIn_;
1677    cost_[sequenceIn_] += modification;
1678    //assert (fabs(modification)<1.0e-6);
1679#ifdef CLP_DEBUG
1680    if ((handler_->logLevel()&32)&&fabs(modification)>1.0e-15)
1681      printf("exact %d new cost %g, change %g\n",sequenceIn_,
1682             cost_[sequenceIn_],modification);
1683#endif
1684#endif
1685   
1686    if (alpha_<0.0) {
1687      // as if from upper bound
1688      directionIn_=-1;
1689      upperIn_=valueIn_;
1690    } else {
1691      // as if from lower bound
1692      directionIn_=1;
1693      lowerIn_=valueIn_;
1694    }
1695  }
1696  if (thisIncrease) {
1697    newTolerance = dualTolerance_+1.0e-4*dblParam_[ClpDualTolerance];
1698  }
1699
1700  // clear arrays
1701
1702  for (i=0;i<2;i++) {
1703    memset(array[i],0,marker[i][0]*sizeof(double));
1704    memset(array[i]+marker[i][1],0,
1705           (numberColumns_-marker[i][1])*sizeof(double));
1706  }
1707}
1708/* Checks if tentative optimal actually means unbounded
1709   Returns -3 if not, 2 if is unbounded */
1710int 
1711ClpSimplexDual::checkUnbounded(CoinIndexedVector * ray,
1712                                   CoinIndexedVector * spare,
1713                                   double changeCost)
1714{
1715  int status=2; // say unbounded
1716  factorization_->updateColumn(spare,ray);
1717  // get reduced cost
1718  int i;
1719  int number=ray->getNumElements();
1720  int * index = ray->getIndices();
1721  double * array = ray->denseVector();
1722  for (i=0;i<number;i++) {
1723    int iRow=index[i];
1724    int iPivot=pivotVariable_[iRow];
1725    changeCost -= cost(iPivot)*array[iRow];
1726  }
1727  double way;
1728  if (changeCost>0.0) {
1729    //try going down
1730    way=1.0;
1731  } else if (changeCost<0.0) {
1732    //try going up
1733    way=-1.0;
1734  } else {
1735#ifdef CLP_DEBUG
1736    printf("can't decide on up or down\n");
1737#endif
1738    way=0.0;
1739    status=-3;
1740  }
1741  double movement=1.0e10*way; // some largish number
1742  double zeroTolerance = 1.0e-14*dualBound_;
1743  for (i=0;i<number;i++) {
1744    int iRow=index[i];
1745    int iPivot=pivotVariable_[iRow];
1746    double arrayValue = array[iRow];
1747    if (fabs(arrayValue)<zeroTolerance)
1748      arrayValue=0.0;
1749    double newValue=solution(iPivot)+movement*arrayValue;
1750    if (newValue>upper(iPivot)+primalTolerance_||
1751        newValue<lower(iPivot)-primalTolerance_)
1752      status=-3; // not unbounded
1753  }
1754  if (status==2) {
1755    // create ray
1756    delete [] ray_;
1757    ray_ = new double [numberColumns_];
1758    ClpFillN(ray_,numberColumns_,0.0);
1759    for (i=0;i<number;i++) {
1760      int iRow=index[i];
1761      int iPivot=pivotVariable_[iRow];
1762      double arrayValue = array[iRow];
1763      if (iPivot<numberColumns_&&fabs(arrayValue)>=zeroTolerance)
1764        ray_[iPivot] = way* array[iRow];
1765    }
1766  }
1767  ray->clear();
1768  return status;
1769}
1770/* Checks if finished.  Updates status */
1771void 
1772ClpSimplexDual::statusOfProblemInDual(int & lastCleaned,int type,
1773                               ClpSimplexProgress &progress)
1774{
1775  if (type==2) {
1776    // trouble - restore solution
1777    memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
1778    memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
1779           numberRows_*sizeof(double));
1780    memcpy(columnActivityWork_,savedSolution_ ,
1781           numberColumns_*sizeof(double));
1782    forceFactorization_=1; // a bit drastic but ..
1783    changeMade_++; // say something changed
1784  }
1785  int tentativeStatus = problemStatus_;
1786  double changeCost;
1787  bool unflagVariables = true;
1788
1789  if (problemStatus_>-3||factorization_->pivots()) {
1790    // factorize
1791    // later on we will need to recover from singularities
1792    // also we could skip if first time
1793    // save dual weights
1794    dualRowPivot_->saveWeights(this,1);
1795    if (type) {
1796      // is factorization okay?
1797      if (internalFactorize(1)) {
1798        // no - restore previous basis
1799        unflagVariables = false;
1800        assert (type==1);
1801        changeMade_++; // say something changed
1802        memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
1803        memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
1804               numberRows_*sizeof(double));
1805        memcpy(columnActivityWork_,savedSolution_ ,
1806               numberColumns_*sizeof(double));
1807        // get correct bounds on all variables
1808        double dummyChangeCost=0.0;
1809        changeBounds(true,rowArray_[2],dummyChangeCost);
1810        // throw away change
1811        int i;
1812        for (i=0;i<4;i++) 
1813          rowArray_[i]->clear();
1814        // need to reject something
1815        char x = isColumn(sequenceOut_) ? 'C' :'R';
1816        handler_->message(CLP_SIMPLEX_FLAG,messages_)
1817          <<x<<sequenceWithin(sequenceOut_)
1818          <<CoinMessageEol;
1819        setFlagged(sequenceOut_);
1820       
1821        forceFactorization_=1; // a bit drastic but ..
1822        type = 2;
1823        //assert (internalFactorize(1)==0);
1824        if (internalFactorize(1)) {
1825          memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
1826          memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
1827                 numberRows_*sizeof(double));
1828          memcpy(columnActivityWork_,savedSolution_ ,
1829                 numberColumns_*sizeof(double));
1830          // debug
1831          assert (internalFactorize(1)==0);
1832        }
1833      }
1834    }
1835    problemStatus_=-3;
1836  }
1837  // at this stage status is -3 or -4 if looks infeasible
1838  // get primal and dual solutions
1839  gutsOfSolution(rowActivityWork_,columnActivityWork_);
1840  // Check if looping
1841  int loop = progress.looping();
1842  int situationChanged=0;
1843  if (loop>=0) {
1844    problemStatus_ = loop; //exit if in loop
1845    if (!problemStatus_) {
1846      // declaring victory
1847      numberPrimalInfeasibilities_ = 0;
1848      sumPrimalInfeasibilities_ = 0.0;
1849    }
1850    return;
1851  } else if (loop<-1) {
1852    // something may have changed
1853    gutsOfSolution(rowActivityWork_,columnActivityWork_);
1854    situationChanged=1;
1855  }
1856  // really for free variables in
1857  if((progressFlag_&2)!=0) {
1858    situationChanged=2;
1859  }
1860  progressFlag_ = 0; //reset progress flag
1861#ifdef CLP_DEBUG
1862  if (!rowScale_&&(handler_->logLevel()&32)) {
1863    double * objectiveSimplex
1864      = ClpCopyOfArray(objective_,numberColumns_,0.0);
1865    double * rowObjectiveSimplex
1866      = ClpCopyOfArray(rowObjective_,numberRows_,0.0);
1867    int i;
1868    double largest;
1869    largest=0.0;
1870    for (i=0;i<numberRows_;i++) {
1871      rowObjectiveSimplex[i] *= optimizationDirection_;
1872      double difference = fabs(rowObjectiveWork_[i]-rowObjectiveSimplex[i]);
1873      if (difference>largest)
1874        largest=difference;
1875    }
1876    for (i=0;i<numberColumns_;i++) {
1877      objectiveSimplex[i] *= optimizationDirection_;
1878      double difference = fabs(objectiveWork_[i]-objectiveSimplex[i]);
1879      if (difference>largest)
1880        largest=difference;
1881    }
1882    if ((handler_->logLevel()&16))
1883      printf("difference in obj %g\n",largest);
1884    delete [] objectiveSimplex;
1885    delete [] rowObjectiveSimplex;
1886  }
1887#endif
1888  handler_->message(CLP_SIMPLEX_STATUS,messages_)
1889    <<numberIterations_<<objectiveValue();
1890  handler_->printing(sumPrimalInfeasibilities_>0.0)
1891    <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
1892  handler_->printing(sumDualInfeasibilities_>0.0)
1893    <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
1894  handler_->printing(numberDualInfeasibilitiesWithoutFree_
1895                     <numberDualInfeasibilities_)
1896                       <<numberDualInfeasibilitiesWithoutFree_;
1897  handler_->message()<<CoinMessageEol;
1898
1899  // dual bound coming in
1900  //double saveDualBound = dualBound_;
1901  while (problemStatus_<=-3) {
1902    bool cleanDuals=situationChanged!=0;
1903    int numberChangedBounds=0;
1904    int doOriginalTolerance=0;
1905    if ( lastCleaned==numberIterations_)
1906      doOriginalTolerance=1;
1907    // check optimal
1908    // give code benefit of doubt
1909    if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
1910        sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
1911      // say optimal (with these bounds etc)
1912      numberDualInfeasibilities_ = 0;
1913      sumDualInfeasibilities_ = 0.0;
1914      numberPrimalInfeasibilities_ = 0;
1915      sumPrimalInfeasibilities_ = 0.0;
1916    }
1917    if (dualFeasible()||problemStatus_==-4) {
1918      if (primalFeasible()) {
1919        // may be optimal - or may be bounds are wrong
1920        handler_->message(CLP_DUAL_BOUNDS,messages_)
1921          <<dualBound_
1922          <<CoinMessageEol;
1923        // save solution in case unbounded
1924        ClpDisjointCopyN(columnActivityWork_,numberColumns_,
1925                          columnArray_[0]->denseVector());
1926        ClpDisjointCopyN(rowActivityWork_,numberRows_,
1927                          rowArray_[2]->denseVector());
1928        numberChangedBounds=changeBounds(false,rowArray_[3],changeCost);
1929        if (numberChangedBounds<=0) {
1930          //looks optimal - do we need to reset tolerance
1931          if (perturbation_==101) {
1932            perturbation_=102; // stop any perturbations
1933            cleanDuals=true;
1934            createRim(4);
1935          }
1936          if (lastCleaned<numberIterations_&&numberTimesOptimal_<4) {
1937            doOriginalTolerance=2;
1938            numberTimesOptimal_++;
1939            changeMade_++; // say something changed
1940            if (numberTimesOptimal_==1) {
1941              dualTolerance_ = min(dualTolerance_,1.0e-8);
1942              // better to have small tolerance even if slower
1943              factorization_->zeroTolerance(1.0e-15);
1944            }
1945            cleanDuals=true;
1946          } else {
1947            problemStatus_=0; // optimal
1948            if (lastCleaned<numberIterations_) {
1949              handler_->message(CLP_SIMPLEX_GIVINGUP,messages_)
1950                <<CoinMessageEol;
1951            }
1952          }
1953        } else {
1954          cleanDuals=true;
1955          if (doOriginalTolerance==1) {
1956            // check unbounded
1957            // find a variable with bad dj
1958            int iSequence;
1959            int iChosen=-1;
1960            double largest = 100.0*primalTolerance_;
1961            for (iSequence=0;iSequence<numberRows_+numberColumns_;
1962                 iSequence++) {
1963              double djValue = dj_[iSequence];
1964              double originalLo = originalLower(iSequence);
1965              double originalUp = originalUpper(iSequence);
1966              if (fabs(djValue)>fabs(largest)) {
1967                if (getStatus(iSequence)!=basic) {
1968                  if (djValue>0&&originalLo<-1.0e20) {
1969                    if (djValue>fabs(largest)) {
1970                      largest=djValue;
1971                      iChosen=iSequence;
1972                    }
1973                  } else if (djValue<0&&originalUp>1.0e20) {
1974                    if (-djValue>fabs(largest)) {
1975                      largest=djValue;
1976                      iChosen=iSequence;
1977                    }
1978                  } 
1979                }
1980              }
1981            }
1982            if (iChosen>=0) {
1983              int iSave=sequenceIn_;
1984              sequenceIn_=iChosen;
1985              unpack(rowArray_[1]);
1986              sequenceIn_ = iSave;
1987              problemStatus_ = checkUnbounded(rowArray_[0],rowArray_[1],
1988                                              changeCost);
1989              rowArray_[1]->clear();
1990            } else {
1991              problemStatus_=-3;
1992            }
1993            if (problemStatus_==2&&perturbation_==101) {
1994              perturbation_=102; // stop any perturbations
1995              cleanDuals=true;
1996              createRim(4);
1997              problemStatus_=-1;
1998            }
1999            if (problemStatus_==2) {
2000              // it is unbounded - restore solution
2001              // but first add in changes to non-basic
2002              int iColumn;
2003              double * original = columnArray_[0]->denseVector();
2004              for (iColumn=0;iColumn<numberColumns_;iColumn++) {
2005                if(getColumnStatus(iColumn)!= basic)
2006                  ray_[iColumn] += 
2007                    columnActivityWork_[iColumn]-original[iColumn];
2008                columnActivityWork_[iColumn] = original[iColumn];
2009              }
2010              ClpDisjointCopyN(rowArray_[2]->denseVector(),numberRows_,
2011                                rowActivityWork_);
2012            }
2013          } else {
2014            doOriginalTolerance=2;
2015            rowArray_[0]->clear();
2016          }
2017        }
2018        ClpFillN(columnArray_[0]->denseVector(),numberColumns_,0.0);
2019        ClpFillN(rowArray_[2]->denseVector(),numberRows_,0.0);
2020      } 
2021      if (problemStatus_==-4) {
2022        // may be infeasible - or may be bounds are wrong
2023        handler_->message(CLP_DUAL_CHECKB,messages_)
2024          <<dualBound_
2025          <<CoinMessageEol;
2026        numberChangedBounds=changeBounds(false,NULL,changeCost);
2027        if (perturbation_==101) {
2028          perturbation_=102; // stop any perturbations
2029          cleanDuals=true;
2030          numberChangedBounds=1;
2031          createRim(4);
2032        }
2033        if (numberChangedBounds<=0||dualBound_>1.0e20||
2034            (largestPrimalError_>1.0&&dualBound_>1.0e17)) {
2035          problemStatus_=1; // infeasible
2036        } else {
2037          problemStatus_=-1; //iterate
2038          cleanDuals=true;
2039          doOriginalTolerance=2;
2040          // and delete ray which has been created
2041          delete [] ray_;
2042          ray_ = NULL;
2043        }
2044
2045      }
2046    } else {
2047      cleanDuals=true;
2048    }
2049    if (problemStatus_<0) {
2050      if (doOriginalTolerance==2) {
2051        // put back original tolerance
2052        lastCleaned=numberIterations_;
2053        handler_->message(CLP_DUAL_ORIGINAL,messages_)
2054          <<CoinMessageEol;
2055
2056        perturbation_=102; // stop any perturbations
2057        createRim(4);
2058        // make sure duals are current
2059        computeDuals();
2060        // put back bounds as they were if was optimal
2061        if (doOriginalTolerance==2) {
2062          changeMade_++; // say something changed
2063          changeBounds(true,NULL,changeCost);
2064          cleanDuals=true;
2065        }
2066      }
2067      if (cleanDuals) {
2068        // make sure dual feasible
2069        // look at all rows and columns
2070        rowArray_[0]->clear();
2071        CoinIotaN(rowArray_[0]->getIndices(),numberRows_,0);
2072        rowArray_[0]->setNumElements(numberRows_);
2073        columnArray_[0]->clear();
2074        CoinIotaN(columnArray_[0]->getIndices(),numberColumns_,0);
2075        columnArray_[0]->setNumElements(numberColumns_);
2076        double objectiveChange=0.0;
2077        updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[1],
2078                          0.0,objectiveChange);
2079        // for now - recompute all
2080        gutsOfSolution(rowActivityWork_,columnActivityWork_);
2081        //assert(numberDualInfeasibilitiesWithoutFree_==0);
2082        if (numberDualInfeasibilities_||situationChanged==2) {
2083          problemStatus_=-1; // carry on as normal
2084        }
2085        situationChanged=0;
2086      } else {
2087        // iterate
2088        problemStatus_=-1;
2089      }
2090    }
2091  }
2092  if (type==0||type==1) {
2093    if (!type) {
2094      // create save arrays
2095      delete [] saveStatus_;
2096      delete [] savedSolution_;
2097      saveStatus_ = new unsigned char [numberRows_+numberColumns_];
2098      savedSolution_ = new double [numberRows_+numberColumns_];
2099    }
2100    // save arrays
2101    memcpy(saveStatus_,status_,(numberColumns_+numberRows_)*sizeof(char));
2102    memcpy(savedSolution_+numberColumns_ ,rowActivityWork_,
2103           numberRows_*sizeof(double));
2104    memcpy(savedSolution_ ,columnActivityWork_,numberColumns_*sizeof(double));
2105  }
2106
2107  // restore weights (if saved) - also recompute infeasibility list
2108  if (tentativeStatus>-3) 
2109    dualRowPivot_->saveWeights(this,(type <2) ? 2 : 4);
2110  else
2111    dualRowPivot_->saveWeights(this,3);
2112  // unflag all variables (we may want to wait a bit?)
2113  if (tentativeStatus!=-2&&unflagVariables) {
2114    int iRow;
2115    for (iRow=0;iRow<numberRows_;iRow++) {
2116      int iPivot=pivotVariable_[iRow];
2117      clearFlagged(iPivot);
2118    }
2119  }
2120  // see if cutoff reached
2121  double limit = 0.0;
2122  getDblParam(ClpDualObjectiveLimit, limit);
2123  if(fabs(limit)<1.0e30&&objectiveValue()*optimizationDirection_>
2124           optimizationDirection_*limit&&
2125           !numberAtFakeBound()) {
2126    problemStatus_=1;
2127  }
2128  if (problemStatus_<0&&!changeMade_) {
2129    problemStatus_=4; // unknown
2130  }
2131  lastGoodIteration_ = numberIterations_;
2132
2133}
2134/* While updateDualsInDual sees what effect is of flip
2135   this does actual flipping.
2136   If change >0.0 then value in array >0.0 => from lower to upper
2137*/
2138void 
2139ClpSimplexDual::flipBounds(CoinIndexedVector * rowArray,
2140                  CoinIndexedVector * columnArray,
2141                  double change)
2142{
2143  double * work;
2144  int number;
2145  int * which;
2146 
2147  int iSection;
2148
2149  for (iSection=0;iSection<2;iSection++) {
2150    int i;
2151    double * solution = solutionRegion(iSection);
2152    double * lower = lowerRegion(iSection);
2153    double * upper = upperRegion(iSection);
2154    int addSequence;
2155    if (!iSection) {
2156      work = rowArray->denseVector();
2157      number = rowArray->getNumElements();
2158      which = rowArray->getIndices();
2159      addSequence = numberColumns_;
2160    } else {
2161      work = columnArray->denseVector();
2162      number = columnArray->getNumElements();
2163      which = columnArray->getIndices();
2164      addSequence = 0;
2165    }
2166   
2167    for (i=0;i<number;i++) {
2168      int iSequence = which[i];
2169#ifndef NDEBUG
2170      double value = work[iSequence]*change;
2171#endif
2172      work[iSequence]=0.0;
2173      Status status = getStatus(iSequence+addSequence);
2174
2175      switch(status) {
2176
2177      case basic:
2178      case isFree:
2179      case superBasic:
2180        break;
2181      case atUpperBound:
2182        assert (value<=0.0);
2183        // to lower bound
2184        setStatus(iSequence+addSequence,atLowerBound);
2185        solution[iSequence] = lower[iSequence];
2186        break;
2187      case atLowerBound:
2188        assert (value>=0.0);
2189        // to upper bound
2190        setStatus(iSequence+addSequence,atUpperBound);
2191        solution[iSequence] = upper[iSequence];
2192        break;
2193      }
2194    }
2195  }
2196  rowArray->setNumElements(0);
2197  columnArray->setNumElements(0);
2198}
2199// Restores bound to original bound
2200void 
2201ClpSimplexDual::originalBound( int iSequence)
2202{
2203  if (iSequence>=numberColumns_) {
2204    // rows
2205    int iRow = iSequence-numberColumns_;
2206    rowLowerWork_[iRow]=rowLower_[iRow];
2207    rowUpperWork_[iRow]=rowUpper_[iRow];
2208    if (rowScale_) {
2209      if (rowLowerWork_[iRow]>-1.0e50)
2210        rowLowerWork_[iRow] *= rowScale_[iRow];
2211      if (rowUpperWork_[iRow]<1.0e50)
2212        rowUpperWork_[iRow] *= rowScale_[iRow];
2213    }
2214  } else {
2215    // columns
2216    columnLowerWork_[iSequence]=columnLower_[iSequence];
2217    columnUpperWork_[iSequence]=columnUpper_[iSequence];
2218    if (rowScale_) {
2219      double multiplier = 1.0/columnScale_[iSequence];
2220      if (columnLowerWork_[iSequence]>-1.0e50)
2221        columnLowerWork_[iSequence] *= multiplier;
2222      if (columnUpperWork_[iSequence]<1.0e50)
2223        columnUpperWork_[iSequence] *= multiplier;
2224    }
2225  }
2226  setFakeBound(iSequence,noFake);
2227}
2228/* As changeBounds but just changes new bounds for a single variable.
2229   Returns true if change */
2230bool 
2231ClpSimplexDual::changeBound( int iSequence)
2232{
2233  // old values
2234  double oldLower=lower_[iSequence];
2235  double oldUpper=upper_[iSequence];
2236  double value=solution_[iSequence];
2237  bool modified=false;
2238  originalBound(iSequence);
2239  // original values
2240  double lowerValue=lower_[iSequence];
2241  double upperValue=upper_[iSequence];
2242  // back to altered values
2243  lower_[iSequence] = oldLower;
2244  upper_[iSequence] = oldUpper;
2245  if (value==oldLower) {
2246    if (upperValue > oldLower + dualBound_) {
2247      upper_[iSequence]=oldLower+dualBound_;
2248      setFakeBound(iSequence,upperFake);
2249      modified=true;
2250    }
2251  } else if (value==oldUpper) {
2252    if (lowerValue < oldUpper - dualBound_) {
2253      lower_[iSequence]=oldUpper-dualBound_;
2254      setFakeBound(iSequence,lowerFake);
2255      modified=true;
2256    }
2257  } else {
2258    assert(value==oldLower||value==oldUpper);
2259  }
2260  return modified;
2261}
2262// Perturbs problem
2263void 
2264ClpSimplexDual::perturb()
2265{
2266  if (perturbation_>100)
2267    return; //perturbed already
2268  int iRow,iColumn;
2269  // dual perturbation
2270  double perturbation=1.0e-20;
2271  // maximum fraction of cost to perturb
2272  double maximumFraction = 1.0e-4;
2273  if (perturbation_>=50) {
2274    perturbation = 1.0e-4;
2275    for (iRow=0;iRow<numberRows_;iRow++) {
2276      double value = fabs(rowObjectiveWork_[iRow]);
2277      perturbation = max(perturbation,value);
2278    }
2279    for (iColumn=0;iColumn<numberColumns_;iColumn++) { 
2280      double value = 
2281        fabs(objectiveWork_[iColumn]);
2282      perturbation = max(perturbation,value);
2283    }
2284    perturbation *= 1.0e-8;
2285  } else if (perturbation_<100) {
2286    perturbation = pow(10.0,perturbation_);
2287    // user is in charge
2288    maximumFraction = 1.0e100;
2289  }
2290  // modify costs
2291  handler_->message(CLP_SIMPLEX_PERTURB,messages_)
2292    <<perturbation
2293    <<CoinMessageEol;
2294  for (iRow=0;iRow<numberRows_;iRow++) {
2295    double value = perturbation;
2296    double currentValue = rowObjectiveWork_[iRow];
2297    value = min(value,maximumFraction*fabs(currentValue)+1.0e-6);
2298    if (rowLowerWork_[iRow]>-largeValue_) {
2299      if (fabs(rowLowerWork_[iRow])<fabs(rowUpperWork_[iRow])) 
2300        value *= CoinDrand48();
2301      else
2302        value *= -CoinDrand48();
2303    } else if (rowUpperWork_[iRow]<largeValue_) {
2304      value *= -CoinDrand48();
2305    } else {
2306      value=0.0;
2307    }
2308    rowObjectiveWork_[iRow] += value;
2309  }
2310  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
2311    double value = perturbation;
2312    double currentValue = objectiveWork_[iColumn];
2313    value = min(value,maximumFraction*fabs(currentValue)+1.0e-6);
2314    if (columnLowerWork_[iColumn]>-largeValue_) {
2315      if (fabs(columnLowerWork_[iColumn])<
2316          fabs(columnUpperWork_[iColumn])) 
2317        value *= CoinDrand48();
2318      else
2319        value *= -CoinDrand48();
2320    } else if (columnUpperWork_[iColumn]<largeValue_) {
2321      value *= -CoinDrand48();
2322    } else {
2323      value=0.0;
2324    }
2325    objectiveWork_[iColumn] += value;
2326  }
2327  // say perturbed
2328  perturbation_=101;
2329
2330}
2331/* For strong branching.  On input lower and upper are new bounds
2332   while on output they are change in objective function values
2333   (>1.0e50 infeasible).
2334   Return code is 0 if nothing interesting, -1 if infeasible both
2335   ways and +1 if infeasible one way (check values to see which one(s))
2336*/
2337int ClpSimplexDual::strongBranching(int numberVariables,const int * variables,
2338                                    double * newLower, double * newUpper,
2339                                    bool stopOnFirstInfeasible,
2340                                    bool alwaysFinish)
2341{
2342  int i;
2343  int returnCode=0;
2344  double saveObjectiveValue = objectiveValue_;
2345#if 1
2346  algorithm_ = -1;
2347
2348  //scaling(false);
2349
2350  // put in standard form (and make row copy)
2351  // create modifiable copies of model rim and do optional scaling
2352  createRim(7+8+16,true);
2353
2354  // change newLower and newUpper if scaled
2355
2356  // Do initial factorization
2357  // and set certain stuff
2358  // We can either set increasing rows so ...IsBasic gives pivot row
2359  // or we can just increment iBasic one by one
2360  // for now let ...iBasic give pivot row
2361  factorization_->increasingRows(2);
2362  // row activities have negative sign
2363  factorization_->slackValue(-1.0);
2364  factorization_->zeroTolerance(1.0e-13);
2365  // save if sparse factorization wanted
2366  int saveSparse = factorization_->sparseThreshold();
2367
2368  int factorizationStatus = internalFactorize(0);
2369  if (factorizationStatus<0)
2370    return 1; // some error
2371  else if (factorizationStatus)
2372    handler_->message(CLP_SINGULARITIES,messages_)
2373      <<factorizationStatus
2374      <<CoinMessageEol;
2375  if (saveSparse) {
2376    // use default at present
2377    factorization_->sparseThreshold(0);
2378    factorization_->goSparse();
2379  }
2380 
2381  // save stuff
2382  ClpFactorization saveFactorization(*factorization_);
2383  // save basis and solution
2384  double * saveSolution = new double[numberRows_+numberColumns_];
2385  memcpy(saveSolution,solution_,
2386         (numberRows_+numberColumns_)*sizeof(double));
2387  unsigned char * saveStatus = 
2388    new unsigned char [numberRows_+numberColumns_];
2389  memcpy(saveStatus,status_,(numberColumns_+numberRows_)*sizeof(char));
2390  // save bounds as createRim makes clean copies
2391  double * saveLower = new double[numberRows_+numberColumns_];
2392  memcpy(saveLower,lower_,
2393         (numberRows_+numberColumns_)*sizeof(double));
2394  double * saveUpper = new double[numberRows_+numberColumns_];
2395  memcpy(saveUpper,upper_,
2396         (numberRows_+numberColumns_)*sizeof(double));
2397  double * saveObjective = new double[numberRows_+numberColumns_];
2398  memcpy(saveObjective,cost_,
2399         (numberRows_+numberColumns_)*sizeof(double));
2400  int * savePivot = new int [numberRows_];
2401  memcpy(savePivot, pivotVariable_, numberRows_*sizeof(int));
2402  // need to save/restore weights.
2403
2404  for (i=0;i<numberVariables;i++) {
2405    int iColumn = variables[i];
2406    double objectiveChange;
2407    double saveBound;
2408   
2409    // try down
2410   
2411    saveBound = columnUpper_[iColumn];
2412    // external view - in case really getting optimal
2413    columnUpper_[iColumn] = newUpper[i];
2414    if (scalingFlag_<=0) 
2415      upper_[iColumn] = newUpper[i];
2416    else 
2417      upper_[iColumn] = newUpper[i]/columnScale_[iColumn]; // scale
2418    // Start of fast iterations
2419    int status = fastDual(alwaysFinish);
2420
2421    // restore
2422    memcpy(solution_,saveSolution,
2423           (numberRows_+numberColumns_)*sizeof(double));
2424    memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char));
2425    memcpy(lower_,saveLower,
2426           (numberRows_+numberColumns_)*sizeof(double));
2427    memcpy(upper_,saveUpper,
2428           (numberRows_+numberColumns_)*sizeof(double));
2429    memcpy(cost_,saveObjective,
2430         (numberRows_+numberColumns_)*sizeof(double));
2431    columnUpper_[iColumn] = saveBound;
2432    memcpy(pivotVariable_, savePivot, numberRows_*sizeof(int));
2433    delete factorization_;
2434    factorization_ = new ClpFactorization(saveFactorization);
2435
2436    if (status||(problemStatus_==0&&!isDualObjectiveLimitReached())) {
2437      objectiveChange = objectiveValue_-saveObjectiveValue;
2438    } else {
2439      objectiveChange = 1.0e100;
2440    }
2441    newUpper[i]=objectiveChange;
2442#ifdef CLP_DEBUG
2443    printf("down on %d costs %g\n",iColumn,objectiveChange);
2444#endif
2445         
2446    // try up
2447   
2448    saveBound = columnLower_[iColumn];
2449    // external view - in case really getting optimal
2450    columnLower_[iColumn] = newLower[i];
2451    if (scalingFlag_<=0) 
2452      lower_[iColumn] = newLower[i];
2453    else 
2454      lower_[iColumn] = newLower[i]/columnScale_[iColumn]; // scale
2455    // Start of fast iterations
2456    status = fastDual(alwaysFinish);
2457
2458    // restore
2459    memcpy(solution_,saveSolution,
2460           (numberRows_+numberColumns_)*sizeof(double));
2461    memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char));
2462    memcpy(lower_,saveLower,
2463           (numberRows_+numberColumns_)*sizeof(double));
2464    memcpy(upper_,saveUpper,
2465           (numberRows_+numberColumns_)*sizeof(double));
2466    memcpy(cost_,saveObjective,
2467         (numberRows_+numberColumns_)*sizeof(double));
2468    columnLower_[iColumn] = saveBound;
2469    memcpy(pivotVariable_, savePivot, numberRows_*sizeof(int));
2470    delete factorization_;
2471    factorization_ = new ClpFactorization(saveFactorization);
2472
2473    if (status||(problemStatus_==0&&!isDualObjectiveLimitReached())) {
2474      objectiveChange = objectiveValue_-saveObjectiveValue;
2475    } else {
2476      objectiveChange = 1.0e100;
2477    }
2478    newLower[i]=objectiveChange;
2479#ifdef CLP_DEBUG
2480    printf("up on %d costs %g\n",iColumn,objectiveChange);
2481#endif
2482         
2483    /* Possibilities are:
2484       Both sides feasible - store
2485       Neither side feasible - set objective high and exit
2486       One side feasible - change bounds and resolve
2487    */
2488    if (newUpper[i]<1.0e100) {
2489      if(newLower[i]<1.0e100) {
2490        // feasible - no action
2491      } else {
2492        // up feasible, down infeasible
2493        returnCode=1;
2494        if (stopOnFirstInfeasible)
2495          break;
2496      }
2497    } else {
2498      if(newLower[i]<1.0e100) {
2499        // down feasible, up infeasible
2500        returnCode=1;
2501        if (stopOnFirstInfeasible)
2502          break;
2503      } else {
2504        // neither side feasible
2505        returnCode=-1;
2506        break;
2507      }
2508    }
2509  }
2510  delete [] saveSolution;
2511  delete [] saveLower;
2512  delete [] saveUpper;
2513  delete [] saveObjective;
2514  delete [] saveStatus;
2515  delete [] savePivot;
2516
2517  factorization_->sparseThreshold(saveSparse);
2518  // Get rid of some arrays and empty factorization
2519  deleteRim();
2520#else
2521  // save basis and solution
2522  double * rowSolution = new double[numberRows_];
2523  memcpy(rowSolution,rowActivity_,numberRows_*sizeof(double));
2524  double * columnSolution = new double[numberColumns_];
2525  memcpy(columnSolution,columnActivity_,numberColumns_*sizeof(double));
2526  unsigned char * saveStatus = 
2527    new unsigned char [numberRows_+numberColumns_];
2528  if (!status_) {
2529    // odd but anyway setup all slack basis
2530    status_ = new unsigned char [numberColumns_+numberRows_];
2531    memset(status_,0,(numberColumns_+numberRows_)*sizeof(char));
2532  }
2533  memcpy(saveStatus,status_,(numberColumns_+numberRows_)*sizeof(char));
2534  for (i=0;i<numberVariables;i++) {
2535    int iColumn = variables[i];
2536    double objectiveChange;
2537   
2538    // try down
2539   
2540    double saveUpper = columnUpper_[iColumn];
2541    columnUpper_[iColumn] = newUpper[i];
2542    dual();
2543
2544    // restore
2545    columnUpper_[iColumn] = saveUpper;
2546    memcpy(rowActivity_,rowSolution,numberRows_*sizeof(double));
2547    memcpy(columnActivity_,columnSolution,numberColumns_*sizeof(double));
2548    memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char));
2549
2550    if (problemStatus_==0&&!isDualObjectiveLimitReached()) {
2551      objectiveChange = objectiveValue_-saveObjectiveValue;
2552    } else {
2553      objectiveChange = 1.0e100;
2554    }
2555    newUpper[i]=objectiveChange;
2556#ifdef CLP_DEBUG
2557    printf("down on %d costs %g\n",iColumn,objectiveChange);
2558#endif
2559         
2560    // try up
2561   
2562    double saveLower = columnLower_[iColumn];
2563    columnLower_[iColumn] = newLower[i];
2564    dual();
2565
2566    // restore
2567    columnLower_[iColumn] = saveLower;
2568    memcpy(rowActivity_,rowSolution,numberRows_*sizeof(double));
2569    memcpy(columnActivity_,columnSolution,numberColumns_*sizeof(double));
2570    memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char));
2571
2572    if (problemStatus_==0&&!isDualObjectiveLimitReached()) {
2573      objectiveChange = objectiveValue_-saveObjectiveValue;
2574    } else {
2575      objectiveChange = 1.0e100;
2576    }
2577    newLower[i]=objectiveChange;
2578#ifdef CLP_DEBUG
2579    printf("up on %d costs %g\n",iColumn,objectiveChange);
2580#endif
2581         
2582    /* Possibilities are:
2583       Both sides feasible - store
2584       Neither side feasible - set objective high and exit
2585       One side feasible - change bounds and resolve
2586    */
2587    if (newUpper[i]<1.0e100) {
2588      if(newLower[i]<1.0e100) {
2589        // feasible - no action
2590      } else {
2591        // up feasible, down infeasible
2592        returnCode=1;
2593        if (stopOnFirstInfeasible)
2594          break;
2595      }
2596    } else {
2597      if(newLower[i]<1.0e100) {
2598        // down feasible, up infeasible
2599        returnCode=1;
2600        if (stopOnFirstInfeasible)
2601          break;
2602      } else {
2603        // neither side feasible
2604        returnCode=-1;
2605        break;
2606      }
2607    }
2608  }
2609  delete [] rowSolution;
2610  delete [] columnSolution;
2611  delete [] saveStatus;
2612#endif
2613  objectiveValue_ = saveObjectiveValue;
2614  return returnCode;
2615}
2616// treat no pivot as finished (unless interesting)
2617int ClpSimplexDual::fastDual(bool alwaysFinish)
2618{
2619  algorithm_ = -1;
2620  dualTolerance_=dblParam_[ClpDualTolerance];
2621  primalTolerance_=dblParam_[ClpPrimalTolerance];
2622
2623  // save dual bound
2624  double saveDualBound = dualBound_;
2625
2626  double objectiveChange;
2627  // for dual we will change bounds using dualBound_
2628  // for this we need clean basis so it is after factorize
2629  gutsOfSolution(rowActivityWork_,columnActivityWork_);
2630  numberFake_ =0; // Number of variables at fake bounds
2631  changeBounds(true,NULL,objectiveChange);
2632
2633  problemStatus_ = -1;
2634  numberIterations_=0;
2635
2636  int lastCleaned=0; // last time objective or bounds cleaned up
2637
2638  // number of times we have declared optimality
2639  numberTimesOptimal_=0;
2640
2641  // Progress indicator
2642  ClpSimplexProgress progress(this);
2643
2644  // This says whether to restore things etc
2645  int factorType=0;
2646  /*
2647    Status of problem:
2648    0 - optimal
2649    1 - infeasible
2650    2 - unbounded
2651    -1 - iterating
2652    -2 - factorization wanted
2653    -3 - redo checking without factorization
2654    -4 - looks infeasible
2655
2656    BUT also from whileIterating return code is:
2657
2658   -1 iterations etc
2659   -2 inaccuracy
2660   -3 slight inaccuracy (and done iterations)
2661   +0 looks optimal (might be unbounded - but we will investigate)
2662   +1 looks infeasible
2663   +3 max iterations
2664
2665  */
2666
2667  int returnCode = 0;
2668
2669  while (problemStatus_<0) {
2670    int iRow,iColumn;
2671    // clear
2672    for (iRow=0;iRow<4;iRow++) {
2673      rowArray_[iRow]->clear();
2674    }   
2675   
2676    for (iColumn=0;iColumn<2;iColumn++) {
2677      columnArray_[iColumn]->clear();
2678    }   
2679
2680    // give matrix (and model costs and bounds a chance to be
2681    // refreshed (normally null)
2682    matrix_->refresh(this);
2683    // may factorize, checks if problem finished
2684    // should be able to speed this up on first time
2685    statusOfProblemInDual(lastCleaned,factorType,progress);
2686
2687    // Say good factorization
2688    factorType=1;
2689
2690    // Do iterations
2691    if (problemStatus_<0) {
2692      returnCode = whileIterating();
2693      if (!alwaysFinish&&returnCode<1) {
2694        double limit = 0.0;
2695        getDblParam(ClpDualObjectiveLimit, limit);
2696        if(fabs(limit)>1.0e30||objectiveValue()*optimizationDirection_<
2697           optimizationDirection_*limit|| 
2698           numberAtFakeBound()) {
2699          // can't say anything interesting - might as well return
2700#ifdef CLP_DEBUG
2701          printf("returning from fastDual after %d iterations with code %d\n",
2702                 numberIterations_,returnCode);
2703#endif
2704          returnCode=1;
2705          break;
2706        }
2707      }
2708      returnCode=0;
2709    }
2710  }
2711
2712  assert(!numberFake_||returnCode||problemStatus_); // all bounds should be okay
2713  dualBound_ = saveDualBound;
2714  return returnCode;
2715}
2716/* Checks number of variables at fake bounds.  This is used by fastDual
2717   so can exit gracefully before end */
2718int 
2719ClpSimplexDual::numberAtFakeBound()
2720{
2721  int iSequence;
2722  int numberFake=0;
2723 
2724  for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
2725    FakeBound bound = getFakeBound(iSequence);
2726    switch(getStatus(iSequence)) {
2727
2728    case basic:
2729      break;
2730    case isFree:
2731    case superBasic:
2732      break;
2733    case atUpperBound:
2734      if (bound==upperFake||bound==bothFake)
2735        numberFake++;
2736      break;
2737    case atLowerBound:
2738      if (bound==lowerFake||bound==bothFake)
2739        numberFake++;
2740      break;
2741    }
2742  }
2743  numberFake_ = numberFake;
2744  return numberFake;
2745}
2746/* Pivot out a variable and choose an incoing one.  Assumes dual
2747   feasible - will not go through a reduced cost. 
2748   Returns step length in theta
2749   Returns ray in ray_ (or NULL if no pivot)
2750   Return codes as before but -1 means no acceptable pivot
2751*/
2752int 
2753ClpSimplexDual::pivotResult()
2754{
2755  abort();
2756  return 0;
2757}
Note: See TracBrowser for help on using the repository browser.