source: tags/release_2002-11-05/ClpSimplexDual.cpp @ 778

Last change on this file since 778 was 50, checked in by ladanyi, 17 years ago

devel-1 merged into HEAD

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