source: branches/devel-1/ClpSimplexDual.cpp @ 2351

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

After more IBM Burlington tests

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 78.1 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        if (whatNext==1) {
621          problemStatus_ =-2; // refactorize
622        } else if (whatNext==2) {
623          // maximum iterations or equivalent
624          problemStatus_= 3;
625          returnCode=3;
626          break;
627        }
628      } else {
629        // no incoming column is valid
630#ifdef CLP_DEBUG
631        if (handler_->logLevel()&32)
632          printf("** no column pivot\n");
633#endif
634        if (factorization_->pivots()<5) {
635          problemStatus_=-4; //say looks infeasible
636          // create ray anyway
637          delete [] ray_;
638          ray_ = new double [ numberRows_];
639          ClpDisjointCopyN(rowArray_[0]->denseVector(),numberRows_,ray_);
640        }
641        rowArray_[0]->clear();
642        columnArray_[0]->clear();
643        returnCode=1;
644        break;
645      }
646    } else {
647      // no pivot row
648#ifdef CLP_DEBUG
649      if (handler_->logLevel()&32)
650        printf("** no row pivot\n");
651#endif
652      if (!factorization_->pivots()) {
653        // may have crept through - so may be optimal
654        //problemStatus_=-5; //say looks unbounded
655        problemStatus_=0;
656        // check any flagged variables
657        int iRow;
658        for (iRow=0;iRow<numberRows_;iRow++) {
659          int iPivot=pivotVariable_[iRow];
660          if (flagged(iPivot))
661            break;
662        }
663        if (iRow<numberRows_) {
664#ifdef CLP_DEBUG
665          std::cerr<<"Flagged variables at end - infeasible?"<<std::endl;
666#endif
667          problemStatus_=-4; //say looks infeasible
668          // create ray anyway
669          delete [] ray_;
670          ray_ = new double [ numberRows_];
671          ClpDisjointCopyN(rowArray_[0]->denseVector(),numberRows_,ray_);
672        }
673      }
674      returnCode=0;
675      break;
676    }
677  }
678  return returnCode;
679}
680/* The duals are updated by the given arrays.
681   Returns number of infeasibilities.
682   rowArray and columnarray will have flipped
683   The output vector has movement (row length array) */
684int
685ClpSimplexDual::updateDualsInDual(CoinIndexedVector * rowArray,
686                                  CoinIndexedVector * columnArray,
687                                  CoinIndexedVector * outputArray,
688                                  double theta,
689                                  double & objectiveChange)
690{
691 
692  outputArray->clear();
693 
694  double * work;
695  int number;
696  int * which;
697 
698  int numberInfeasibilities=0;
699  int numberRowInfeasibilities=0;
700 
701  // see whether we will be doing full recompute
702  bool fullRecompute= (rowArray->getNumElements()==numberRows_&&
703                       columnArray->getNumElements()==numberColumns_);
704  int numberAtFake=0;
705 
706  // use a tighter tolerance except for all being okay
707  double tolerance = dualTolerance_;
708 
709  double changeObj=0.0;
710 
711  int iSection;
712 
713  for (iSection=0;iSection<2;iSection++) {
714    int i;
715    double * solution = solutionRegion(iSection);
716    double * reducedCost = djRegion(iSection);
717    double * lower = lowerRegion(iSection);
718    double * upper = upperRegion(iSection);
719    double * cost = costRegion(iSection);
720    int addSequence;
721    if (!iSection) {
722      addSequence = numberColumns_;
723      work = rowArray->denseVector();
724      number = rowArray->getNumElements();
725      which = rowArray->getIndices();
726    } else {
727      // set number of infeasibilities in row array
728      addSequence=0;
729      numberRowInfeasibilities=numberInfeasibilities;
730      rowArray->setNumElements(numberInfeasibilities);
731      numberInfeasibilities=0;
732      work = columnArray->denseVector();
733      number = columnArray->getNumElements();
734      which = columnArray->getIndices();
735    }
736   
737    for (i=0;i<number;i++) {
738      int iSequence = which[i];
739      double alphaI = work[iSequence];
740      double value = reducedCost[iSequence]-theta*alphaI;
741      work[iSequence]=0.0;
742      reducedCost[iSequence]=value;
743     
744      if (!fixed(iSequence+addSequence)) {
745        double movement=0.0;
746        FakeBound bound = getFakeBound(iSequence+addSequence);
747        Status status = getStatus(iSequence+addSequence);
748
749        switch(status) {
750         
751        case basic:
752        case superBasic:
753          break;
754        case isFree:
755          if (fabs(value)>tolerance) { 
756#ifdef CLP_DEBUG
757            if (handler_->logLevel()&32)
758              printf("%d %d, free has dj of %g, alpha %g\n",
759                     iSection,iSequence,value,alphaI);
760#endif
761          }
762          break;
763        case atUpperBound:
764          if (value>tolerance) {
765            // to lower bound (if swap)
766            // put back alpha
767            which[numberInfeasibilities++]=iSequence;
768            work[iSequence]=alphaI;
769            movement = lower[iSequence]-upper[iSequence];
770#ifdef CLP_DEBUG
771            if ((handler_->logLevel()&32))
772              printf("%d %d, new dj %g, alpha %g, movement %g\n",
773                   iSection,iSequence,value,alphaI,movement);
774#endif
775            changeObj += movement*cost[iSequence];
776            if (bound==ClpSimplexDual::bothFake||
777                bound==ClpSimplexDual::lowerFake) 
778              numberAtFake++;
779          } else if (fullRecompute) {
780            // at correct bound
781            if (bound==ClpSimplexDual::bothFake||
782                bound==ClpSimplexDual::upperFake) {
783              // but flip if dj would allow
784              if (bound==ClpSimplexDual::upperFake&&
785                  value>=-tolerance) {
786                movement = lower[iSequence]-upper[iSequence];
787                setStatus(iSequence+addSequence,atLowerBound);
788                solution[iSequence] = lower[iSequence];
789                changeObj += movement*cost[iSequence];
790              } else {
791                numberAtFake++;
792              }
793            }
794          }
795          break;
796        case atLowerBound:
797          if (value<-tolerance) {
798            // to upper bound
799            // put back alpha
800            which[numberInfeasibilities++]=iSequence;
801            work[iSequence]=alphaI;
802            movement = upper[iSequence] - lower[iSequence];
803#ifdef CLP_DEBUG
804            if ((handler_->logLevel()&32))
805              printf("%d %d, new dj %g, alpha %g, movement %g\n",
806                   iSection,iSequence,value,alphaI,movement);
807#endif
808            changeObj += movement*cost[iSequence];
809            if (bound==ClpSimplexDual::bothFake||
810                bound==ClpSimplexDual::upperFake) 
811              numberAtFake++;
812          } else if (fullRecompute) {
813            // at correct bound
814            if (bound==ClpSimplexDual::bothFake||
815                bound==ClpSimplexDual::lowerFake) {
816              // but flip if dj would allow
817              if (bound==ClpSimplexDual::lowerFake&&
818                  value<=tolerance) {
819                movement = upper[iSequence] - lower[iSequence];
820                setStatus(iSequence+addSequence,atUpperBound);
821                solution[iSequence] = upper[iSequence];
822                changeObj += movement*cost[iSequence];
823              } else {
824                numberAtFake++;
825              }
826            }
827          }
828          break;
829        }
830        if (!fullRecompute) {
831          if (movement) {
832            if (!iSection) {
833              // row (sign ?)
834              outputArray->quickAdd(iSequence,-movement);
835            } else {
836              matrix_->add(this,outputArray,iSequence,movement);
837            }
838          }
839        }
840      }
841    }
842  }
843#ifdef CLP_DEBUG
844  if (fullRecompute&&numberAtFake&&(handler_->logLevel()&16)!=0) 
845    printf("%d fake after full update\n",numberAtFake);
846#endif
847  outputArray->stopQuickAdd();
848  // set number of infeasibilities
849  columnArray->setNumElements(numberInfeasibilities);
850  numberInfeasibilities += numberRowInfeasibilities;
851  if (fullRecompute) {
852    // do actual flips
853    flipBounds(rowArray,columnArray,theta);
854    numberFake_ = numberAtFake;
855  }
856  objectiveChange += changeObj;
857  return numberInfeasibilities;
858}
859/*
860   Chooses dual pivot row
861   Would be faster with separate region to scan
862   and will have this (with square of infeasibility) when steepest
863   For easy problems we can just choose one of the first rows we look at
864*/
865void
866ClpSimplexDual::dualRow()
867{
868  // get pivot row using whichever method it is
869  pivotRow_=dualRowPivot_->pivotRow();
870  if (pivotRow_>=0) {
871    int iPivot=pivotVariable_[pivotRow_];
872    sequenceOut_ = iPivot;
873    if (iPivot>=numberColumns_) {
874      // slack
875      iPivot-=numberColumns_;
876      valueOut_=rowActivityWork_[iPivot];
877      lowerOut_=rowLowerWork_[iPivot];
878      upperOut_=rowUpperWork_[iPivot];
879    } else {
880      // column
881      valueOut_=columnActivityWork_[iPivot];
882      lowerOut_=columnLowerWork_[iPivot];
883      upperOut_=columnUpperWork_[iPivot];
884    }
885    // if we have problems we could try other way and hope we get a
886    // zero pivot?
887    if (valueOut_>upperOut_) {
888      directionOut_ = -1;
889      dualOut_ = valueOut_ - upperOut_;
890    } else {
891      directionOut_ = 1;
892      dualOut_ = lowerOut_ - valueOut_;
893    }
894#ifdef CLP_DEBUG
895    assert(dualOut_>=0.0);
896#endif
897  }
898  return ;
899}
900// Checks if any fake bounds active - if so returns number and modifies
901// dualBound_ and everything.
902// Free variables will be left as free
903// Returns number of bounds changed if >=0
904// Returns -1 if not initialize and no effect
905// Fills in changeVector which can be used to see if unbounded
906// and cost of change vector
907int
908ClpSimplexDual::changeBounds(bool initialize,
909                                 CoinIndexedVector * outputArray,
910                                 double & changeCost)
911{
912  if (!initialize) {
913    int numberInfeasibilities;
914    double newBound;
915    newBound = 5.0*dualBound_;
916    numberInfeasibilities=0;
917    changeCost=0.0;
918    // put back original bounds and then check
919    createRim(3);
920    int iSequence;
921    // bounds will get bigger - just look at ones at bounds
922    for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
923      double lowerValue=lower_[iSequence];
924      double upperValue=upper_[iSequence];
925      double value=solution_[iSequence];
926      setFakeBound(iSequence,ClpSimplexDual::noFake);
927      switch(getStatus(iSequence)) {
928
929      case basic:
930        break;
931      case isFree:
932      case superBasic:
933        break;
934      case atUpperBound:
935        if (fabs(value-upperValue)>primalTolerance_) 
936          numberInfeasibilities++;
937        break;
938      case atLowerBound:
939        if (fabs(value-lowerValue)>primalTolerance_) 
940          numberInfeasibilities++;
941        break;
942      }
943    }
944    if (numberInfeasibilities) {
945      int iSequence;
946      for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
947        double lowerValue=lower_[iSequence];
948        double upperValue=upper_[iSequence];
949        double newLowerValue;
950        double newUpperValue;
951        Status status = getStatus(iSequence);
952        if (status==atUpperBound||
953            status==atLowerBound) {
954          double value = solution_[iSequence];
955          if (value-lowerValue<=upperValue-value) {
956            newLowerValue = max(lowerValue,value-0.666667*newBound);
957            newUpperValue = min(upperValue,newLowerValue+newBound);
958          } else {
959            newUpperValue = min(upperValue,value+0.666667*newBound);
960            newLowerValue = max(lowerValue,newUpperValue-newBound);
961          }
962          lower_[iSequence]=newLowerValue;
963          upper_[iSequence]=newUpperValue;
964          if (newLowerValue > lowerValue) {
965            if (newUpperValue < upperValue) 
966              setFakeBound(iSequence,ClpSimplexDual::bothFake);
967            else
968              setFakeBound(iSequence,ClpSimplexDual::lowerFake);
969          } else {
970            if (newUpperValue < upperValue) 
971              setFakeBound(iSequence,ClpSimplexDual::upperFake);
972          }
973          if (status==atUpperBound)
974            solution_[iSequence] = newUpperValue;
975          else
976            solution_[iSequence] = newLowerValue;
977          double movement = solution_[iSequence] - value;
978          if (movement&&outputArray) {
979            if (iSequence>=numberColumns_) {
980              outputArray->quickAdd(iSequence,-movement);
981              changeCost += movement*cost_[iSequence];
982            } else {
983              matrix_->add(this,outputArray,iSequence,movement);
984              changeCost += movement*cost_[iSequence];
985            }
986          }
987        }
988      }
989      dualBound_ = newBound;
990      if (outputArray)
991        outputArray->stopQuickAdd();
992    } else {
993      numberInfeasibilities=-1;
994    }
995    return numberInfeasibilities;
996  } else {
997    int iSequence;
998     
999    for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
1000      Status status = getStatus(iSequence);
1001      if (status==atUpperBound||
1002          status==atLowerBound) {
1003        double lowerValue=lower_[iSequence];
1004        double upperValue=upper_[iSequence];
1005        double value = solution_[iSequence];
1006        if (lowerValue>-largeValue_||upperValue<largeValue_) {
1007          if (lowerValue-value>-0.5*dualBound_||
1008              upperValue-value<0.5*dualBound_) {
1009            if (fabs(lowerValue-value)<=fabs(upperValue-value)) {
1010              if (upperValue > lowerValue + dualBound_) {
1011                upper_[iSequence]=lowerValue+dualBound_;
1012                setFakeBound(iSequence,ClpSimplexDual::upperFake);
1013              }
1014            } else {
1015              if (lowerValue < upperValue - dualBound_) {
1016                lower_[iSequence]=upperValue-dualBound_;
1017                setFakeBound(iSequence,ClpSimplexDual::lowerFake);
1018              }
1019            }
1020          } else {
1021            lower_[iSequence]=-0.5*dualBound_;
1022            upper_[iSequence]= 0.5*dualBound_;
1023            setFakeBound(iSequence,ClpSimplexDual::bothFake);
1024          }
1025        }
1026      }
1027    }
1028    return 1;
1029  }
1030}
1031/*
1032   Row array has row part of pivot row (as duals so sign may be switched)
1033   Column array has column part.
1034   This chooses pivot column.
1035   Spare array will be needed when we start getting clever.
1036   We will check for basic so spare array will never overflow.
1037   If necessary will modify costs
1038*/
1039void
1040ClpSimplexDual::dualColumn(CoinIndexedVector * rowArray,
1041                       CoinIndexedVector * columnArray,
1042                       CoinIndexedVector * spareArray,
1043                       CoinIndexedVector * spareArray2)
1044{
1045  double * work;
1046  int number;
1047  int * which;
1048  double * reducedCost;
1049
1050  int iSection;
1051
1052  sequenceIn_=-1;
1053  int numberPossiblySwapped=0;
1054  int numberRemaining=0;
1055 
1056  double totalThru=0.0; // for when variables flip
1057  double acceptablePivot=1.0e-7;
1058  if (factorization_->pivots())
1059    acceptablePivot=1.0e-5; // if we have iterated be more strict
1060  double bestEverPivot=acceptablePivot;
1061  int lastSequence = -1;
1062  double lastPivot=0.0;
1063  double upperTheta;
1064  double newTolerance = dualTolerance_;
1065  // will we need to increase tolerance
1066  bool thisIncrease=false;
1067  // If we think we need to modify costs (not if something from broad sweep)
1068  bool modifyCosts=false;
1069  // Increase in objective due to swapping bounds (may be negative)
1070  double increaseInObjective=0.0;
1071
1072  // use spareArrays to put ones looked at in
1073  // we are going to flip flop between
1074  int iFlip = 0;
1075  // Possible list of pivots
1076  int interesting[2];
1077  // where possible swapped ones are
1078  int swapped[2];
1079  // for zeroing out arrays after
1080  int marker[2][2];
1081  // pivot elements
1082  double * array[2], * spare, * spare2;
1083  // indices
1084  int * indices[2], * index, * index2;
1085  spareArray->clear();
1086  spareArray2->clear();
1087  array[0] = spareArray->denseVector();
1088  indices[0] = spareArray->getIndices();
1089  spare = array[0];
1090  index = indices[0];
1091  array[1] = spareArray2->denseVector();
1092  indices[1] = spareArray2->getIndices();
1093  int i;
1094  double * lower;
1095  double * upper;
1096
1097  // initialize lists
1098  for (i=0;i<2;i++) {
1099    interesting[i]=0;
1100    swapped[i]=numberColumns_;
1101    marker[i][0]=0;
1102    marker[i][1]=numberColumns_;
1103  }
1104
1105  /*
1106    First we get a list of possible pivots.  We can also see if the
1107    problem looks infeasible or whether we want to pivot in free variable.
1108    This may make objective go backwards but can only happen a finite
1109    number of times and I do want free variables basic.
1110
1111    Then we flip back and forth.  At the start of each iteration
1112    interesting[iFlip] should have possible candidates and swapped[iFlip]
1113    will have pivots if we decide to take a previous pivot.
1114    At end of each iteration interesting[1-iFlip] should have
1115    candidates if we go through this theta and swapped[1-iFlip]
1116    pivots if we don't go through.
1117
1118    At first we increase theta and see what happens.  We start
1119    theta at a reasonable guess.  If in right area then we do bit by bit.
1120
1121   */
1122
1123  // do first pass to get possibles
1124  // We can also see if infeasible or pivoting on free
1125  double tentativeTheta = 1.0e22;
1126  upperTheta = 1.0e31;
1127  double freePivot = acceptablePivot;
1128  for (iSection=0;iSection<2;iSection++) {
1129
1130    int addSequence;
1131
1132    if (!iSection) {
1133      lower = rowLowerWork_;
1134      upper = rowUpperWork_;
1135      work = rowArray->denseVector();
1136      number = rowArray->getNumElements();
1137      which = rowArray->getIndices();
1138      reducedCost = rowReducedCost_;
1139      addSequence = numberColumns_;
1140    } else {
1141      lower = columnLowerWork_;
1142      upper = columnUpperWork_;
1143      work = columnArray->denseVector();
1144      number = columnArray->getNumElements();
1145      which = columnArray->getIndices();
1146      reducedCost = reducedCostWork_;
1147      addSequence = 0;
1148    }
1149   
1150    for (i=0;i<number;i++) {
1151      int iSequence = which[i];
1152      double alpha = work[iSequence];
1153      if (fixed(iSequence+addSequence)||!alpha) 
1154        continue; // skip fixed ones or (zeroed out)
1155      double oldValue = reducedCost[iSequence];
1156      double value = oldValue-tentativeTheta*alpha;
1157      int keep = 0;
1158
1159      switch(getStatus(iSequence+addSequence)) {
1160         
1161      case basic:
1162        break;
1163      case isFree:
1164      case superBasic:
1165        if (oldValue>dualTolerance_) {
1166          if (value<-newTolerance) 
1167            keep = 2;
1168        } else if (oldValue<-dualTolerance_) {
1169          if (value>newTolerance) 
1170            keep = 2;
1171        } else {
1172          if (alpha>=acceptablePivot) 
1173            keep = 2;
1174          else if (-alpha>=acceptablePivot) 
1175            keep = 2;
1176        }
1177        break;
1178      case atUpperBound:
1179#ifdef CHECK_DJ
1180        // For Debug so we can find out where problem is
1181        perturbation_ = iSequence+addSequence;
1182#endif
1183        assert (oldValue<=dualTolerance_*1.0001);
1184        if (value>newTolerance) {
1185          keep = 1;
1186          value = oldValue-upperTheta*alpha;
1187          if (value>newTolerance && -alpha>=acceptablePivot) 
1188            upperTheta = (oldValue-newTolerance)/alpha;
1189        }
1190        break;
1191      case atLowerBound:
1192#ifdef CHECK_DJ
1193        // For Debug so we can find out where problem is
1194        perturbation_ = iSequence+addSequence;
1195#endif
1196        assert (oldValue>=-dualTolerance_*1.0001);
1197        if (value<-newTolerance) {
1198          keep = 1;
1199          value = oldValue-upperTheta*alpha;
1200          if (value<-newTolerance && alpha>=acceptablePivot) 
1201            upperTheta = (oldValue+newTolerance)/alpha;
1202        }
1203        break;
1204      }
1205      if (keep) {
1206        if (keep==2) {
1207          // free - choose largest
1208          if (fabs(alpha)>freePivot) {
1209            freePivot=fabs(alpha);
1210            sequenceIn_ = iSequence + addSequence;
1211            theta_=oldValue/alpha;
1212          }
1213        } else {
1214          // add to list
1215          spare[numberRemaining]=alpha;
1216          index[numberRemaining++]=iSequence+addSequence;
1217        }
1218      }
1219    }
1220  }
1221  interesting[0]=numberRemaining;
1222  marker[0][0] = numberRemaining;
1223
1224  if (!numberRemaining)
1225    return; // Looks infeasible
1226
1227  if (sequenceIn_>=0) {
1228    // free variable - always choose
1229  } else {
1230
1231    theta_=1.0e50;
1232    // now flip flop between spare arrays until reasonable theta
1233    tentativeTheta = max(10.0*upperTheta,1.0e-7);
1234   
1235    // loops increasing tentative theta until can't go through
1236   
1237    while (tentativeTheta < 1.0e22) {
1238      double thruThis = 0.0;
1239     
1240      double bestPivot=acceptablePivot;
1241      int bestSequence=-1;
1242     
1243      numberPossiblySwapped = numberColumns_;
1244      numberRemaining = 0;
1245
1246      upperTheta = 1.0e50;
1247
1248      spare = array[iFlip];
1249      index = indices[iFlip];
1250      spare2 = array[1-iFlip];
1251      index2 = indices[1-iFlip];
1252
1253      // try 3 different ways
1254      // 1 bias increase by ones with slightly wrong djs
1255      // 2 bias by all
1256      // 3 bias by all - tolerance (doesn't seem very good)
1257#define TRYBIAS 1
1258
1259
1260      double increaseInThis=0.0; //objective increase in this loop
1261     
1262      for (i=0;i<interesting[iFlip];i++) {
1263        int iSequence = index[i];
1264        double alpha = spare[i];
1265        double oldValue = dj_[iSequence];
1266        double value = oldValue-tentativeTheta*alpha;
1267
1268        if (alpha < 0.0) {
1269          //at upper bound
1270          if (value>newTolerance) {
1271            double range = upper_[iSequence] - lower_[iSequence];
1272            thruThis -= range*alpha;
1273#if TRYBIAS==1
1274            if (oldValue>0.0)
1275              increaseInThis -= oldValue*range;
1276#elif TRYBIAS==2
1277            increaseInThis -= oldValue*range;
1278#else
1279            increaseInThis -= (oldValue+dualTolerance_)*range;
1280#endif
1281            // goes on swapped list (also means candidates if too many)
1282            spare2[--numberPossiblySwapped]=alpha;
1283            index2[numberPossiblySwapped]=iSequence;
1284            if (fabs(alpha)>bestPivot) {
1285              bestPivot=fabs(alpha);
1286              bestSequence=numberPossiblySwapped;
1287            }
1288          } else {
1289            value = oldValue-upperTheta*alpha;
1290            if (value>newTolerance && -alpha>=acceptablePivot) 
1291              upperTheta = (oldValue-newTolerance)/alpha;
1292            spare2[numberRemaining]=alpha;
1293            index2[numberRemaining++]=iSequence;
1294          }
1295        } else {
1296          // at lower bound
1297          if (value<-newTolerance) {
1298            double range = upper_[iSequence] - lower_[iSequence];
1299            thruThis += range*alpha;
1300            //?? is this correct - and should we look at good ones
1301#if TRYBIAS==1
1302            if (oldValue<0.0)
1303              increaseInThis += oldValue*range;
1304#elif TRYBIAS==2
1305            increaseInThis += oldValue*range;
1306#else
1307            increaseInThis += (oldValue-dualTolerance_)*range;
1308#endif
1309            // goes on swapped list (also means candidates if too many)
1310            spare2[--numberPossiblySwapped]=alpha;
1311            index2[numberPossiblySwapped]=iSequence;
1312            if (fabs(alpha)>bestPivot) {
1313              bestPivot=fabs(alpha);
1314              bestSequence=numberPossiblySwapped;
1315            }
1316          } else {
1317            value = oldValue-upperTheta*alpha;
1318            if (value<-newTolerance && alpha>=acceptablePivot) 
1319              upperTheta = (oldValue+newTolerance)/alpha;
1320            spare2[numberRemaining]=alpha;
1321            index2[numberRemaining++]=iSequence;
1322          }
1323        }
1324      }
1325      swapped[1-iFlip]=numberPossiblySwapped;
1326      interesting[1-iFlip]=numberRemaining;
1327      marker[1-iFlip][0]= max(marker[1-iFlip][0],numberRemaining);
1328      marker[1-iFlip][1]= min(marker[1-iFlip][1],numberPossiblySwapped);
1329     
1330      if (totalThru+thruThis>=fabs(dualOut_)||
1331          increaseInObjective+increaseInThis<0.0) {
1332        // We should be pivoting in this batch
1333        // so compress down to this lot
1334        numberRemaining=0;
1335        for (i=numberColumns_-1;i>=swapped[1-iFlip];i--) {
1336          spare[numberRemaining]=spare2[i];
1337          index[numberRemaining++]=index2[i];
1338        }
1339        interesting[iFlip]=numberRemaining;
1340        int iTry;
1341#define MAXTRY 100
1342        // first get ratio with tolerance
1343        for (iTry=0;iTry<MAXTRY;iTry++) {
1344         
1345          upperTheta=1.0e50;
1346          numberPossiblySwapped = numberColumns_;
1347          numberRemaining = 0;
1348
1349          increaseInThis=0.0; //objective increase in this loop
1350
1351          thruThis=0.0;
1352
1353          spare = array[iFlip];
1354          index = indices[iFlip];
1355          spare2 = array[1-iFlip];
1356          index2 = indices[1-iFlip];
1357     
1358          for (i=0;i<interesting[iFlip];i++) {
1359            int iSequence=index[i];
1360            double alpha=spare[i];
1361            double oldValue = dj_[iSequence];
1362            double value = oldValue-upperTheta*alpha;
1363           
1364            if (alpha < 0.0) {
1365              //at upper bound
1366              if (value>newTolerance) {
1367                if (-alpha>=acceptablePivot) {
1368                  upperTheta = (oldValue-newTolerance)/alpha;
1369                }
1370              }
1371            } else {
1372              // at lower bound
1373              if (value<-newTolerance) {
1374                if (alpha>=acceptablePivot) {
1375                  upperTheta = (oldValue+newTolerance)/alpha;
1376                }
1377              }
1378            }
1379          }
1380          bestPivot=acceptablePivot;
1381          sequenceIn_=-1;
1382          // now choose largest and sum all ones which will go through
1383#define MINIMUMTHETA 1.0e-12
1384          for (i=0;i<interesting[iFlip];i++) {
1385            int iSequence=index[i];
1386            double alpha=spare[i];
1387            double value = dj_[iSequence]-upperTheta*alpha;
1388            double badDj=0.0;
1389           
1390            bool addToSwapped=false;
1391           
1392            if (alpha < 0.0) {
1393              //at upper bound
1394              if (value>=0.0) { 
1395                addToSwapped=true;
1396#if TRYBIAS==1
1397                badDj = -max(dj_[iSequence],0.0);
1398#elif TRYBIAS==2
1399                badDj = -dj_[iSequence];
1400#else
1401                badDj = -dj_[iSequence]-dualTolerance_;
1402#endif
1403              }
1404            } else {
1405              // at lower bound
1406              if (value<=0.0) {
1407                addToSwapped=true;
1408#if TRYBIAS==1
1409                badDj = min(dj_[iSequence],0.0);
1410#elif TRYBIAS==2
1411                badDj = dj_[iSequence];
1412#else
1413                badDj = dj_[iSequence]-dualTolerance_;
1414#endif
1415              }
1416            }
1417            if (!addToSwapped) {
1418              // add to list of remaining
1419              spare2[numberRemaining]=alpha;
1420              index2[numberRemaining++]=iSequence;
1421            } else {
1422              // add to list of swapped
1423              spare2[--numberPossiblySwapped]=alpha;
1424              index2[numberPossiblySwapped]=iSequence;
1425              // select if largest pivot
1426              if (fabs(alpha)>bestPivot) {
1427                sequenceIn_ = numberPossiblySwapped;
1428                bestPivot =  fabs(alpha);
1429                theta_ = dj_[iSequence]/alpha;
1430              }
1431              double range = upper[iSequence] - lower[iSequence];
1432              thruThis += range*fabs(alpha);
1433              increaseInThis += badDj*range;
1434            }
1435          }
1436          swapped[1-iFlip]=numberPossiblySwapped;
1437          interesting[1-iFlip]=numberRemaining;
1438          marker[1-iFlip][0]= max(marker[1-iFlip][0],numberRemaining);
1439          marker[1-iFlip][1]= min(marker[1-iFlip][1],numberPossiblySwapped);
1440          // If we stop now this will be increase in objective (I think)
1441          double increase = (fabs(dualOut_)-totalThru)*theta_;
1442          increase += increaseInObjective;
1443          if (theta_<0.0)
1444            thruThis += fabs(dualOut_); // force using this one
1445          if (increaseInObjective<0.0&&increase<0.0&&lastSequence>=0) {
1446            // back
1447            // We may need to be more careful - we could do by
1448            // switch so we always do fine grained?
1449            bestPivot=0.0;
1450          } else {
1451            // add in
1452            totalThru += thruThis;
1453            increaseInObjective += increaseInThis;
1454          }
1455          if (bestPivot<0.1*bestEverPivot&&
1456              bestEverPivot>1.0e-6&&bestPivot<1.0e-3) {
1457            // back to previous one
1458            sequenceIn_=lastSequence;
1459            // swap regions
1460            iFlip = 1-iFlip;
1461            break;
1462          } else if (sequenceIn_==-1&&upperTheta>largeValue_) {
1463            if (lastPivot>acceptablePivot) {
1464              // back to previous one
1465              sequenceIn_=lastSequence;
1466              // swap regions
1467              iFlip = 1-iFlip;
1468            } else {
1469              // can only get here if all pivots too small
1470            }
1471            break;
1472          } else if (totalThru>=fabs(dualOut_)) {
1473            modifyCosts=true; // fine grain - we can modify costs
1474            break; // no point trying another loop
1475          } else {
1476            lastSequence=sequenceIn_;
1477            if (bestPivot>bestEverPivot)
1478              bestEverPivot=bestPivot;
1479            iFlip = 1 -iFlip;
1480            modifyCosts=true; // fine grain - we can modify costs
1481        }
1482        }
1483        if (iTry==MAXTRY)
1484          iFlip = 1-iFlip; // flip back
1485        break;
1486      } else {
1487        // skip this lot
1488        if (bestPivot>1.0e-3||bestPivot>bestEverPivot) {
1489          bestEverPivot=bestPivot;
1490          lastSequence=bestSequence;
1491        } else {
1492          // keep old swapped
1493          memcpy(array[1-iFlip]+swapped[iFlip],
1494                 array[iFlip]+swapped[iFlip],
1495                 (numberColumns_-swapped[iFlip])*sizeof(double));
1496          memcpy(indices[1-iFlip]+swapped[iFlip],
1497                 indices[iFlip]+swapped[iFlip],
1498                 (numberColumns_-swapped[iFlip])*sizeof(int));
1499          marker[1-iFlip][1] = min(marker[1-iFlip][1],swapped[iFlip]);
1500          swapped[1-iFlip]=swapped[iFlip];
1501        }
1502        increaseInObjective += increaseInThis;
1503        iFlip = 1 - iFlip; // swap regions
1504        tentativeTheta = 2.0*upperTheta;
1505        totalThru += thruThis;
1506      }
1507    }
1508   
1509    // can get here without sequenceIn_ set but with lastSequence
1510    if (sequenceIn_<0&&lastSequence>=0) {
1511      // back to previous one
1512      sequenceIn_=lastSequence;
1513      // swap regions
1514      iFlip = 1-iFlip;
1515    }
1516   
1517    if (sequenceIn_>=0) {
1518      // at this stage sequenceIn_ is just pointer into index array
1519      // flip just so we can use iFlip
1520      iFlip = 1 -iFlip;
1521      spare = array[iFlip];
1522      index = indices[iFlip];
1523      double oldValue;
1524      double alpha = spare[sequenceIn_];
1525      sequenceIn_ = indices[iFlip][sequenceIn_];
1526      oldValue = dj_[sequenceIn_];
1527      theta_ = oldValue/alpha;
1528#if 0
1529      if (numberIterations_>2000)
1530        exit(99);
1531      if (numberIterations_>2000-20)
1532        handler_->setLogLevel(63);
1533      if (numberIterations_>2000-20)
1534        printf("theta %g oldValue %g tol %g %g\n",theta_,oldValue,dualTolerance_,
1535               newTolerance);
1536#endif
1537      if (theta_<MINIMUMTHETA) {
1538        // can't pivot to zero
1539        if (oldValue-MINIMUMTHETA*alpha>=-dualTolerance_) {
1540          theta_=MINIMUMTHETA;
1541        } else if (oldValue-MINIMUMTHETA*alpha>=-newTolerance) {
1542          theta_=MINIMUMTHETA;
1543          thisIncrease=true;
1544        } else {
1545          theta_=(oldValue+newTolerance)/alpha;
1546          assert(theta_>=0.0);
1547          thisIncrease=true;
1548        } 
1549      }
1550      // may need to adjust costs so all dual feasible AND pivoted is exactly 0
1551      if (modifyCosts) {
1552        int i;
1553        double * workRow = rowArray->denseVector();
1554        double * workColumn = columnArray->denseVector();
1555        for (i=numberColumns_-1;i>=swapped[iFlip];i--) {
1556          int iSequence=index[i];
1557          double alpha;
1558          if (iSequence>=numberColumns_)
1559            alpha=workRow[iSequence-numberColumns_];
1560          else
1561            alpha=workColumn[iSequence];
1562          double value = dj_[iSequence]-theta_*alpha;
1563#if 0
1564          if (numberIterations_>2000-20)
1565            printf("%d alpha %g value %g\n",iSequence,alpha,value);
1566#endif
1567           
1568          // can't be free here
1569         
1570          if (alpha < 0.0) {
1571            //at upper bound
1572            if (value>dualTolerance_) {
1573              thisIncrease=true;
1574#define MODIFYCOST 2
1575#if MODIFYCOST
1576              // modify cost to hit new tolerance
1577              double modification = alpha*theta_-dj_[iSequence]
1578                +newTolerance;
1579              //modification = min(modification,dualTolerance_);
1580              //assert (fabs(modification)<1.0e-7);
1581              dj_[iSequence] += modification;
1582              cost_[iSequence] +=  modification;
1583#if 0
1584              if (numberIterations_>2000-20)
1585                printf("%d acost %g mod %g\n",iSequence,cost_[iSequence],
1586                       modification);
1587#endif
1588#endif
1589            }
1590          } else {
1591            // at lower bound
1592            if (-value>dualTolerance_) {
1593              thisIncrease=true;
1594#if MODIFYCOST
1595              // modify cost to hit new tolerance
1596              double modification = alpha*theta_-dj_[iSequence]
1597                -newTolerance;
1598              //modification = max(modification,-dualTolerance_);
1599              //assert (fabs(modification)<1.0e-7);
1600              dj_[iSequence] += modification;
1601              cost_[iSequence] +=  modification;
1602#if 0
1603              if (numberIterations_>2000-20)
1604                printf("%d bcost %g mod %g\n",iSequence,cost_[iSequence],
1605                       modification);
1606#endif
1607#endif
1608            }
1609          }
1610        }
1611      }
1612    }
1613  }
1614
1615  if (sequenceIn_>=0) {
1616    if (sequenceIn_>=numberColumns_) {
1617      //slack
1618      alpha_ = rowArray->denseVector()[sequenceIn_-numberColumns_];
1619    } else {
1620      // column
1621      alpha_ = columnArray->denseVector()[sequenceIn_];
1622    }
1623    lowerIn_ = lower_[sequenceIn_];
1624    upperIn_ = upper_[sequenceIn_];
1625    valueIn_ = solution_[sequenceIn_];
1626    dualIn_ = dj_[sequenceIn_];
1627
1628    if (numberTimesOptimal_) {
1629      // can we adjust cost back closer to original
1630      //*** add coding
1631    }
1632#if MODIFYCOST>1   
1633    // modify cost to hit zero exactly
1634    // so (dualIn_+modification)==theta_*alpha_
1635    double modification = theta_*alpha_-dualIn_;
1636    dualIn_ += modification;
1637    dj_[sequenceIn_]=dualIn_;
1638    cost_[sequenceIn_] += modification;
1639    //assert (fabs(modification)<1.0e-6);
1640#ifdef CLP_DEBUG
1641    if ((handler_->logLevel()&32)&&fabs(modification)>1.0e-15)
1642      printf("exact %d new cost %g, change %g\n",sequenceIn_,
1643             cost_[sequenceIn_],modification);
1644#endif
1645#endif
1646   
1647    if (alpha_<0.0) {
1648      // as if from upper bound
1649      directionIn_=-1;
1650      upperIn_=valueIn_;
1651    } else {
1652      // as if from lower bound
1653      directionIn_=1;
1654      lowerIn_=valueIn_;
1655    }
1656  }
1657  if (thisIncrease) {
1658    newTolerance = dualTolerance_+1.0e-4*dblParam_[ClpDualTolerance];
1659  }
1660
1661  // clear arrays
1662
1663  for (i=0;i<2;i++) {
1664    memset(array[i],0,marker[i][0]*sizeof(double));
1665    memset(array[i]+marker[i][1],0,
1666           (numberColumns_-marker[i][1])*sizeof(double));
1667  }
1668}
1669/* Checks if tentative optimal actually means unbounded
1670   Returns -3 if not, 2 if is unbounded */
1671int 
1672ClpSimplexDual::checkUnbounded(CoinIndexedVector * ray,
1673                                   CoinIndexedVector * spare,
1674                                   double changeCost)
1675{
1676  int status=2; // say unbounded
1677  factorization_->updateColumn(spare,ray);
1678  // get reduced cost
1679  int i;
1680  int number=ray->getNumElements();
1681  int * index = ray->getIndices();
1682  double * array = ray->denseVector();
1683  for (i=0;i<number;i++) {
1684    int iRow=index[i];
1685    int iPivot=pivotVariable_[iRow];
1686    changeCost -= cost(iPivot)*array[iRow];
1687  }
1688  double way;
1689  if (changeCost>0.0) {
1690    //try going down
1691    way=1.0;
1692  } else if (changeCost<0.0) {
1693    //try going up
1694    way=-1.0;
1695  } else {
1696#ifdef CLP_DEBUG
1697    printf("can't decide on up or down\n");
1698#endif
1699    way=0.0;
1700    status=-3;
1701  }
1702  double movement=1.0e10*way; // some largish number
1703  double zeroTolerance = 1.0e-14*dualBound_;
1704  for (i=0;i<number;i++) {
1705    int iRow=index[i];
1706    int iPivot=pivotVariable_[iRow];
1707    double arrayValue = array[iRow];
1708    if (fabs(arrayValue)<zeroTolerance)
1709      arrayValue=0.0;
1710    double newValue=solution(iPivot)+movement*arrayValue;
1711    if (newValue>upper(iPivot)+primalTolerance_||
1712        newValue<lower(iPivot)-primalTolerance_)
1713      status=-3; // not unbounded
1714  }
1715  if (status==2) {
1716    // create ray
1717    delete [] ray_;
1718    ray_ = new double [numberColumns_];
1719    ClpFillN(ray_,numberColumns_,0.0);
1720    for (i=0;i<number;i++) {
1721      int iRow=index[i];
1722      int iPivot=pivotVariable_[iRow];
1723      double arrayValue = array[iRow];
1724      if (iPivot<numberColumns_&&fabs(arrayValue)>=zeroTolerance)
1725        ray_[iPivot] = way* array[iRow];
1726    }
1727  }
1728  ray->clear();
1729  return status;
1730}
1731/* Checks if finished.  Updates status */
1732void 
1733ClpSimplexDual::statusOfProblemInDual(int & lastCleaned,int type,
1734                               ClpSimplexProgress &progress)
1735{
1736  if (type==2) {
1737    // trouble - restore solution
1738    memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
1739    memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
1740           numberRows_*sizeof(double));
1741    memcpy(columnActivityWork_,savedSolution_ ,
1742           numberColumns_*sizeof(double));
1743    forceFactorization_=1; // a bit drastic but ..
1744    changeMade_++; // say something changed
1745  }
1746  int tentativeStatus = problemStatus_;
1747  double changeCost;
1748
1749  if (problemStatus_>-3) {
1750    // factorize
1751    // later on we will need to recover from singularities
1752    // also we could skip if first time
1753    // save dual weights
1754    dualRowPivot_->saveWeights(this,1);
1755    if (type) {
1756      // is factorization okay?
1757      if (internalFactorize(1)) {
1758        // no - restore previous basis
1759        assert (type==1);
1760        changeMade_++; // say something changed
1761        memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
1762        memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
1763               numberRows_*sizeof(double));
1764        memcpy(columnActivityWork_,savedSolution_ ,
1765               numberColumns_*sizeof(double));
1766        // get correct bounds on all variables
1767        double dummyChangeCost=0.0;
1768        changeBounds(true,rowArray_[2],dummyChangeCost);
1769        // throw away change
1770        rowArray_[2]->clear();
1771        forceFactorization_=1; // a bit drastic but ..
1772        type = 2;
1773        assert (internalFactorize(1)==0);
1774      }
1775    }
1776    problemStatus_=-3;
1777  }
1778  // at this stage status is -3 or -4 if looks infeasible
1779  // get primal and dual solutions
1780  gutsOfSolution(rowActivityWork_,columnActivityWork_);
1781  // Check if looping
1782  int loop = progress.looping();
1783  bool situationChanged=false;
1784  if (loop>=0) {
1785    problemStatus_ = loop; //exit if in loop
1786    return;
1787  } else if (loop<-1) {
1788    // something may have changed
1789    gutsOfSolution(rowActivityWork_,columnActivityWork_);
1790    situationChanged=true;
1791  }
1792  progressFlag_ = 0; //reset progress flag
1793#ifdef CLP_DEBUG
1794  if (!rowScale_&&(handler_->logLevel()&32)) {
1795    double * objectiveSimplex
1796      = ClpCopyOfArray(objective_,numberColumns_,0.0);
1797    double * rowObjectiveSimplex
1798      = ClpCopyOfArray(rowObjective_,numberRows_,0.0);
1799    int i;
1800    double largest;
1801    largest=0.0;
1802    for (i=0;i<numberRows_;i++) {
1803      rowObjectiveSimplex[i] *= optimizationDirection_;
1804      double difference = fabs(rowObjectiveWork_[i]-rowObjectiveSimplex[i]);
1805      if (difference>largest)
1806        largest=difference;
1807    }
1808    for (i=0;i<numberColumns_;i++) {
1809      objectiveSimplex[i] *= optimizationDirection_;
1810      double difference = fabs(objectiveWork_[i]-objectiveSimplex[i]);
1811      if (difference>largest)
1812        largest=difference;
1813    }
1814    if ((handler_->logLevel()&16))
1815      printf("difference in obj %g\n",largest);
1816    delete [] objectiveSimplex;
1817    delete [] rowObjectiveSimplex;
1818  }
1819#endif
1820  handler_->message(CLP_SIMPLEX_STATUS,messages_)
1821    <<numberIterations_<<objectiveValue();
1822  handler_->printing(sumPrimalInfeasibilities_>0.0)
1823    <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
1824  handler_->printing(sumDualInfeasibilities_>0.0)
1825    <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
1826  handler_->printing(numberDualInfeasibilitiesWithoutFree_
1827                     <numberDualInfeasibilities_)
1828                       <<numberDualInfeasibilities_-
1829    numberDualInfeasibilitiesWithoutFree_;
1830  handler_->message()<<CoinMessageEol;
1831
1832  while (problemStatus_<=-3) {
1833    bool cleanDuals=situationChanged;
1834    int numberChangedBounds=0;
1835    int doOriginalTolerance=0;
1836    if ( lastCleaned==numberIterations_)
1837      doOriginalTolerance=1;
1838    // check optimal
1839    // give code benefit of doubt
1840    if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
1841        sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
1842      // say optimal (with these bounds etc)
1843      numberDualInfeasibilities_ = 0;
1844      sumDualInfeasibilities_ = 0.0;
1845      numberPrimalInfeasibilities_ = 0;
1846      sumPrimalInfeasibilities_ = 0.0;
1847    }
1848    if (dualFeasible()||problemStatus_==-4) {
1849      if (primalFeasible()) {
1850        // may be optimal - or may be bounds are wrong
1851        handler_->message(CLP_DUAL_BOUNDS,messages_)
1852          <<dualBound_
1853          <<CoinMessageEol;
1854        // save solution in case unbounded
1855        ClpDisjointCopyN(columnActivityWork_,numberColumns_,
1856                          columnArray_[0]->denseVector());
1857        ClpDisjointCopyN(rowActivityWork_,numberRows_,
1858                          rowArray_[2]->denseVector());
1859        numberChangedBounds=changeBounds(false,rowArray_[0],changeCost);
1860        if (numberChangedBounds<=0) {
1861          //looks optimal - do we need to reset tolerance
1862          if (lastCleaned<numberIterations_&&numberTimesOptimal_<4) {
1863            doOriginalTolerance=2;
1864            numberTimesOptimal_++;
1865            changeMade_++; // say something changed
1866            if (numberTimesOptimal_==1) {
1867              dualTolerance_ = min(dualTolerance_,1.0e-8);
1868              // better to have small tolerance even if slower
1869              factorization_->zeroTolerance(1.0e-15);
1870            }
1871            cleanDuals=true;
1872          } else {
1873            problemStatus_=0; // optimal
1874            if (lastCleaned<numberIterations_) {
1875              handler_->message(CLP_SIMPLEX_GIVINGUP,messages_)
1876                <<CoinMessageEol;
1877            }
1878          }
1879        } else {
1880          cleanDuals=true;
1881          if (doOriginalTolerance==1) {
1882            // check unbounded
1883            problemStatus_ = checkUnbounded(rowArray_[0],rowArray_[1],
1884                                            changeCost);
1885            if (problemStatus_==2) {
1886              // it is unbounded - restore solution
1887              // but first add in changes to non-basic
1888              int iColumn;
1889              double * original = columnArray_[0]->denseVector();
1890              for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1891                if(getColumnStatus(iColumn)!= basic)
1892                  ray_[iColumn] += 
1893                    columnActivityWork_[iColumn]-original[iColumn];
1894                columnActivityWork_[iColumn] = original[iColumn];
1895              }
1896              ClpDisjointCopyN(rowArray_[2]->denseVector(),numberRows_,
1897                                rowActivityWork_);
1898            }
1899          } else {
1900            doOriginalTolerance=2;
1901            rowArray_[0]->clear();
1902          }
1903        }
1904        ClpFillN(columnArray_[0]->denseVector(),numberColumns_,0.0);
1905        ClpFillN(rowArray_[2]->denseVector(),numberRows_,0.0);
1906      } 
1907      if (problemStatus_==-4) {
1908        // may be infeasible - or may be bounds are wrong
1909        handler_->message(CLP_DUAL_CHECKB,messages_)
1910          <<dualBound_
1911          <<CoinMessageEol;
1912        numberChangedBounds=changeBounds(false,NULL,changeCost);
1913        if (numberChangedBounds<=0||dualBound_>1.0e20||
1914            (largestPrimalError_>1.0&&dualBound_>1.0e17)) {
1915          problemStatus_=1; // infeasible
1916        } else {
1917          problemStatus_=-1; //iterate
1918          cleanDuals=true;
1919          doOriginalTolerance=2;
1920          // and delete ray which has been created
1921          delete [] ray_;
1922          ray_ = NULL;
1923        }
1924
1925      }
1926    } else {
1927      cleanDuals=true;
1928    }
1929    if (problemStatus_<0) {
1930      if (doOriginalTolerance==2) {
1931        // put back original tolerance
1932        lastCleaned=numberIterations_;
1933        handler_->message(CLP_DUAL_ORIGINAL,messages_)
1934          <<CoinMessageEol;
1935
1936        perturbation_=102; // stop any perturbations
1937        createRim(4);
1938        // make sure duals are current
1939        computeDuals();
1940        // put back bounds as they were if was optimal
1941        if (doOriginalTolerance==2) {
1942          changeMade_++; // say something changed
1943          changeBounds(true,NULL,changeCost);
1944          cleanDuals=true;
1945        }
1946      }
1947      if (cleanDuals) {
1948        // make sure dual feasible
1949        // look at all rows and columns
1950        rowArray_[0]->clear();
1951        CoinIotaN(rowArray_[0]->getIndices(),numberRows_,0);
1952        rowArray_[0]->setNumElements(numberRows_);
1953        columnArray_[0]->clear();
1954        CoinIotaN(columnArray_[0]->getIndices(),numberColumns_,0);
1955        columnArray_[0]->setNumElements(numberColumns_);
1956        double objectiveChange=0.0;
1957        updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[1],
1958                          0.0,objectiveChange);
1959        // for now - recompute all
1960        gutsOfSolution(rowActivityWork_,columnActivityWork_);
1961        assert(numberDualInfeasibilitiesWithoutFree_==0);
1962        if (numberDualInfeasibilities_) {
1963          // bad free variables
1964          if (primalFeasible()) {
1965            std::cerr<<"Free variable problem?"<<std::endl;
1966            abort(); // what now
1967          }
1968          problemStatus_=-1; // carry on as normal
1969        }
1970      } else {
1971        // iterate
1972        problemStatus_=-1;
1973      }
1974    }
1975  }
1976  if (type==0||type==1) {
1977    if (!type) {
1978      // create save arrays
1979      delete [] saveStatus_;
1980      delete [] savedSolution_;
1981      saveStatus_ = new unsigned char [numberRows_+numberColumns_];
1982      savedSolution_ = new double [numberRows_+numberColumns_];
1983    }
1984    // save arrays
1985    memcpy(saveStatus_,status_,(numberColumns_+numberRows_)*sizeof(char));
1986    memcpy(savedSolution_+numberColumns_ ,rowActivityWork_,
1987           numberRows_*sizeof(double));
1988    memcpy(savedSolution_ ,columnActivityWork_,numberColumns_*sizeof(double));
1989  }
1990
1991  // restore weights (if saved) - also recompute infeasibility list
1992  if (tentativeStatus>-3) 
1993    dualRowPivot_->saveWeights(this,(type <2) ? 2 : 4);
1994  else
1995    dualRowPivot_->saveWeights(this,3);
1996  // unflag all variables (we may want to wait a bit?)
1997  int iRow;
1998  for (iRow=0;iRow<numberRows_;iRow++) {
1999    int iPivot=pivotVariable_[iRow];
2000    clearFlagged(iPivot);
2001  }
2002  if (problemStatus_<0&&!changeMade_) {
2003    problemStatus_=4; // unknown
2004  }
2005
2006}
2007/* While updateDualsInDual sees what effect is of flip
2008   this does actual flipping.
2009   If change >0.0 then value in array >0.0 => from lower to upper
2010*/
2011void 
2012ClpSimplexDual::flipBounds(CoinIndexedVector * rowArray,
2013                  CoinIndexedVector * columnArray,
2014                  double change)
2015{
2016  double * work;
2017  int number;
2018  int * which;
2019 
2020  int iSection;
2021
2022  for (iSection=0;iSection<2;iSection++) {
2023    int i;
2024    double * solution = solutionRegion(iSection);
2025    double * lower = lowerRegion(iSection);
2026    double * upper = upperRegion(iSection);
2027    int addSequence;
2028    if (!iSection) {
2029      work = rowArray->denseVector();
2030      number = rowArray->getNumElements();
2031      which = rowArray->getIndices();
2032      addSequence = numberColumns_;
2033    } else {
2034      work = columnArray->denseVector();
2035      number = columnArray->getNumElements();
2036      which = columnArray->getIndices();
2037      addSequence = 0;
2038    }
2039   
2040    for (i=0;i<number;i++) {
2041      int iSequence = which[i];
2042#ifndef NDEBUG
2043      double value = work[iSequence]*change;
2044#endif
2045      work[iSequence]=0.0;
2046      Status status = getStatus(iSequence+addSequence);
2047
2048      switch(status) {
2049
2050      case basic:
2051      case isFree:
2052      case superBasic:
2053        break;
2054      case atUpperBound:
2055        assert (value<=0.0);
2056        // to lower bound
2057        setStatus(iSequence+addSequence,atLowerBound);
2058        solution[iSequence] = lower[iSequence];
2059        break;
2060      case atLowerBound:
2061        assert (value>=0.0);
2062        // to upper bound
2063        setStatus(iSequence+addSequence,atUpperBound);
2064        solution[iSequence] = upper[iSequence];
2065        break;
2066      }
2067    }
2068  }
2069  rowArray->setNumElements(0);
2070  columnArray->setNumElements(0);
2071}
2072// Restores bound to original bound
2073void 
2074ClpSimplexDual::originalBound( int iSequence)
2075{
2076  if (iSequence>=numberColumns_) {
2077    // rows
2078    int iRow = iSequence-numberColumns_;
2079    rowLowerWork_[iRow]=rowLower_[iRow];
2080    rowUpperWork_[iRow]=rowUpper_[iRow];
2081    if (rowScale_) {
2082      if (rowLowerWork_[iRow]>-1.0e50)
2083        rowLowerWork_[iRow] *= rowScale_[iRow];
2084      if (rowUpperWork_[iRow]<1.0e50)
2085        rowUpperWork_[iRow] *= rowScale_[iRow];
2086    }
2087  } else {
2088    // columns
2089    columnLowerWork_[iSequence]=columnLower_[iSequence];
2090    columnUpperWork_[iSequence]=columnUpper_[iSequence];
2091    if (rowScale_) {
2092      double multiplier = 1.0/columnScale_[iSequence];
2093      if (columnLowerWork_[iSequence]>-1.0e50)
2094        columnLowerWork_[iSequence] *= multiplier;
2095      if (columnUpperWork_[iSequence]<1.0e50)
2096        columnUpperWork_[iSequence] *= multiplier;
2097    }
2098  }
2099  setFakeBound(iSequence,noFake);
2100}
2101/* As changeBounds but just changes new bounds for a single variable.
2102   Returns true if change */
2103bool 
2104ClpSimplexDual::changeBound( int iSequence)
2105{
2106  // old values
2107  double oldLower=lower_[iSequence];
2108  double oldUpper=upper_[iSequence];
2109  double value=solution_[iSequence];
2110  bool modified=false;
2111  originalBound(iSequence);
2112  // original values
2113  double lowerValue=lower_[iSequence];
2114  double upperValue=upper_[iSequence];
2115  // back to altered values
2116  lower_[iSequence] = oldLower;
2117  upper_[iSequence] = oldUpper;
2118  if (value==oldLower) {
2119    if (upperValue > oldLower + dualBound_) {
2120      upper_[iSequence]=oldLower+dualBound_;
2121      setFakeBound(iSequence,upperFake);
2122      modified=true;
2123    }
2124  } else if (value==oldUpper) {
2125    if (lowerValue < oldUpper - dualBound_) {
2126      lower_[iSequence]=oldUpper-dualBound_;
2127      setFakeBound(iSequence,lowerFake);
2128      modified=true;
2129    }
2130  } else {
2131    assert(value==oldLower||value==oldUpper);
2132  }
2133  return modified;
2134}
2135// Perturbs problem
2136void 
2137ClpSimplexDual::perturb()
2138{
2139  if (perturbation_>100)
2140    return; //perturbed already
2141  int iRow,iColumn;
2142  // dual perturbation
2143  double perturbation=1.0e-20;
2144  // maximum fraction of cost to perturb
2145  double maximumFraction = 1.0e-4;
2146  if (perturbation_==100) {
2147    perturbation = 1.0e-4;
2148    for (iRow=0;iRow<numberRows_;iRow++) {
2149      double value = fabs(rowActivityWork_[iRow]*rowObjectiveWork_[iRow]);
2150      perturbation = max(perturbation,value);
2151    }
2152    for (iColumn=0;iColumn<numberColumns_;iColumn++) { 
2153      double value = 
2154        fabs(columnActivityWork_[iColumn]*objectiveWork_[iColumn]);
2155      perturbation = max(perturbation,value);
2156    }
2157    perturbation *= 1.0e-8;
2158  } else if (perturbation_<100) {
2159    perturbation = pow(10.0,perturbation_);
2160    // user is in charge
2161    maximumFraction = 1.0e100;
2162  }
2163  // modify costs
2164  handler_->message(CLP_SIMPLEX_PERTURB,messages_)
2165    <<perturbation
2166    <<CoinMessageEol;
2167  for (iRow=0;iRow<numberRows_;iRow++) {
2168    double value = perturbation;
2169    double currentValue = rowObjectiveWork_[iRow];
2170    value = min(value,maximumFraction*fabs(currentValue)+1.0e-6);
2171    if (rowLowerWork_[iRow]>-largeValue_) {
2172      if (fabs(rowLowerWork_[iRow])<fabs(rowUpperWork_[iRow])) 
2173        value *= drand48();
2174      else
2175        value *= -drand48();
2176    } else if (rowUpperWork_[iRow]<largeValue_) {
2177      value *= -drand48();
2178    } else {
2179      value=0.0;
2180    }
2181    rowObjectiveWork_[iRow] += value;
2182  }
2183  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
2184    double value = perturbation;
2185    double currentValue = objectiveWork_[iColumn];
2186    value = min(value,maximumFraction*fabs(currentValue)+1.0e-6);
2187    if (columnLowerWork_[iColumn]>-largeValue_) {
2188      if (fabs(columnLowerWork_[iColumn])<
2189          fabs(columnUpperWork_[iColumn])) 
2190        value *= drand48();
2191      else
2192        value *= -drand48();
2193    } else if (columnUpperWork_[iColumn]<largeValue_) {
2194      value *= -drand48();
2195    } else {
2196      value=0.0;
2197    }
2198    objectiveWork_[iColumn] += value;
2199  }
2200  // say perturbed
2201  perturbation_=102;
2202
2203}
2204/* For strong branching.  On input lower and upper are new bounds
2205   while on output they are change in objective function values
2206   (>1.0e50 infeasible).
2207   Return code is 0 if nothing interesting, -1 if infeasible both
2208   ways and +1 if infeasible one way (check values to see which one(s))
2209*/
2210int ClpSimplexDual::strongBranching(int numberVariables,const int * variables,
2211                                    double * newLower, double * newUpper,
2212                                    bool stopOnFirstInfeasible,
2213                                    bool alwaysFinish)
2214{
2215  int i;
2216  int returnCode=0;
2217  double saveObjectiveValue = objectiveValue_;
2218#if 1
2219  algorithm_ = -1;
2220
2221  //scaling(false);
2222
2223  // put in standard form (and make row copy)
2224  // create modifiable copies of model rim and do optional scaling
2225  createRim(7+8+16,true);
2226
2227  // change newLower and newUpper if scaled
2228
2229  // Do initial factorization
2230  // and set certain stuff
2231  // We can either set increasing rows so ...IsBasic gives pivot row
2232  // or we can just increment iBasic one by one
2233  // for now let ...iBasic give pivot row
2234  factorization_->increasingRows(2);
2235  // row activities have negative sign
2236  factorization_->slackValue(-1.0);
2237  factorization_->zeroTolerance(1.0e-13);
2238  // save if sparse factorization wanted
2239  int saveSparse = factorization_->sparseThreshold();
2240
2241  int factorizationStatus = internalFactorize(0);
2242  if (factorizationStatus<0)
2243    return 1; // some error
2244  else if (factorizationStatus)
2245    handler_->message(CLP_SINGULARITIES,messages_)
2246      <<factorizationStatus
2247      <<CoinMessageEol;
2248  if (saveSparse) {
2249    // use default at present
2250    factorization_->sparseThreshold(0);
2251    factorization_->goSparse();
2252  }
2253 
2254  // save stuff
2255  ClpFactorization saveFactorization(*factorization_);
2256  // save basis and solution
2257  double * saveSolution = new double[numberRows_+numberColumns_];
2258  memcpy(saveSolution,solution_,
2259         (numberRows_+numberColumns_)*sizeof(double));
2260  unsigned char * saveStatus = 
2261    new unsigned char [numberRows_+numberColumns_];
2262  memcpy(saveStatus,status_,(numberColumns_+numberRows_)*sizeof(char));
2263  // save bounds as createRim makes clean copies
2264  double * saveLower = new double[numberRows_+numberColumns_];
2265  memcpy(saveLower,lower_,
2266         (numberRows_+numberColumns_)*sizeof(double));
2267  double * saveUpper = new double[numberRows_+numberColumns_];
2268  memcpy(saveUpper,upper_,
2269         (numberRows_+numberColumns_)*sizeof(double));
2270  double * saveObjective = new double[numberRows_+numberColumns_];
2271  memcpy(saveObjective,cost_,
2272         (numberRows_+numberColumns_)*sizeof(double));
2273  int * savePivot = new int [numberRows_];
2274  memcpy(savePivot, pivotVariable_, numberRows_*sizeof(int));
2275  // need to save/restore weights.
2276
2277  for (i=0;i<numberVariables;i++) {
2278    int iColumn = variables[i];
2279    double objectiveChange;
2280    double saveBound;
2281   
2282    // try down
2283   
2284    saveBound = columnUpper_[iColumn];
2285    // external view - in case really getting optimal
2286    columnUpper_[iColumn] = newUpper[i];
2287    if (scalingFlag_<=0) 
2288      upper_[iColumn] = newUpper[i];
2289    else 
2290      upper_[iColumn] = newUpper[i]/columnScale_[iColumn]; // scale
2291    // Start of fast iterations
2292    int status = fastDual(alwaysFinish);
2293
2294    // restore
2295    memcpy(solution_,saveSolution,
2296           (numberRows_+numberColumns_)*sizeof(double));
2297    memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char));
2298    memcpy(lower_,saveLower,
2299           (numberRows_+numberColumns_)*sizeof(double));
2300    memcpy(upper_,saveUpper,
2301           (numberRows_+numberColumns_)*sizeof(double));
2302    memcpy(cost_,saveObjective,
2303         (numberRows_+numberColumns_)*sizeof(double));
2304    columnUpper_[iColumn] = saveBound;
2305    memcpy(pivotVariable_, savePivot, numberRows_*sizeof(int));
2306    delete factorization_;
2307    factorization_ = new ClpFactorization(saveFactorization);
2308
2309    if (status||(problemStatus_==0&&!isDualObjectiveLimitReached())) {
2310      objectiveChange = objectiveValue_-saveObjectiveValue;
2311    } else {
2312      objectiveChange = 1.0e100;
2313    }
2314    newUpper[i]=objectiveChange;
2315#ifdef CLP_DEBUG
2316    printf("down on %d costs %g\n",iColumn,objectiveChange);
2317#endif
2318         
2319    // try up
2320   
2321    saveBound = columnLower_[iColumn];
2322    // external view - in case really getting optimal
2323    columnLower_[iColumn] = newLower[i];
2324    if (scalingFlag_<=0) 
2325      lower_[iColumn] = newLower[i];
2326    else 
2327      lower_[iColumn] = newLower[i]/columnScale_[iColumn]; // scale
2328    // Start of fast iterations
2329    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    columnLower_[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    newLower[i]=objectiveChange;
2352#ifdef CLP_DEBUG
2353    printf("up on %d costs %g\n",iColumn,objectiveChange);
2354#endif
2355         
2356    /* Possibilities are:
2357       Both sides feasible - store
2358       Neither side feasible - set objective high and exit
2359       One side feasible - change bounds and resolve
2360    */
2361    if (newUpper[i]<1.0e100) {
2362      if(newLower[i]<1.0e100) {
2363        // feasible - no action
2364      } else {
2365        // up feasible, down infeasible
2366        returnCode=1;
2367        if (stopOnFirstInfeasible)
2368          break;
2369      }
2370    } else {
2371      if(newLower[i]<1.0e100) {
2372        // down feasible, up infeasible
2373        returnCode=1;
2374        if (stopOnFirstInfeasible)
2375          break;
2376      } else {
2377        // neither side feasible
2378        returnCode=-1;
2379        break;
2380      }
2381    }
2382  }
2383  delete [] saveSolution;
2384  delete [] saveLower;
2385  delete [] saveUpper;
2386  delete [] saveObjective;
2387  delete [] saveStatus;
2388  delete [] savePivot;
2389
2390  factorization_->sparseThreshold(saveSparse);
2391  // Get rid of some arrays and empty factorization
2392  deleteRim();
2393#else
2394  // save basis and solution
2395  double * rowSolution = new double[numberRows_];
2396  memcpy(rowSolution,rowActivity_,numberRows_*sizeof(double));
2397  double * columnSolution = new double[numberColumns_];
2398  memcpy(columnSolution,columnActivity_,numberColumns_*sizeof(double));
2399  unsigned char * saveStatus = 
2400    new unsigned char [numberRows_+numberColumns_];
2401  if (!status_) {
2402    // odd but anyway setup all slack basis
2403    status_ = new unsigned char [numberColumns_+numberRows_];
2404    memset(status_,0,(numberColumns_+numberRows_)*sizeof(char));
2405  }
2406  memcpy(saveStatus,status_,(numberColumns_+numberRows_)*sizeof(char));
2407  for (i=0;i<numberVariables;i++) {
2408    int iColumn = variables[i];
2409    double objectiveChange;
2410   
2411    // try down
2412   
2413    double saveUpper = columnUpper_[iColumn];
2414    columnUpper_[iColumn] = newUpper[i];
2415    dual();
2416
2417    // restore
2418    columnUpper_[iColumn] = saveUpper;
2419    memcpy(rowActivity_,rowSolution,numberRows_*sizeof(double));
2420    memcpy(columnActivity_,columnSolution,numberColumns_*sizeof(double));
2421    memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char));
2422
2423    if (problemStatus_==0&&!isDualObjectiveLimitReached()) {
2424      objectiveChange = objectiveValue_-saveObjectiveValue;
2425    } else {
2426      objectiveChange = 1.0e100;
2427    }
2428    newUpper[i]=objectiveChange;
2429#ifdef CLP_DEBUG
2430    printf("down on %d costs %g\n",iColumn,objectiveChange);
2431#endif
2432         
2433    // try up
2434   
2435    double saveLower = columnLower_[iColumn];
2436    columnLower_[iColumn] = newLower[i];
2437    dual();
2438
2439    // restore
2440    columnLower_[iColumn] = saveLower;
2441    memcpy(rowActivity_,rowSolution,numberRows_*sizeof(double));
2442    memcpy(columnActivity_,columnSolution,numberColumns_*sizeof(double));
2443    memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char));
2444
2445    if (problemStatus_==0&&!isDualObjectiveLimitReached()) {
2446      objectiveChange = objectiveValue_-saveObjectiveValue;
2447    } else {
2448      objectiveChange = 1.0e100;
2449    }
2450    newLower[i]=objectiveChange;
2451#ifdef CLP_DEBUG
2452    printf("up on %d costs %g\n",iColumn,objectiveChange);
2453#endif
2454         
2455    /* Possibilities are:
2456       Both sides feasible - store
2457       Neither side feasible - set objective high and exit
2458       One side feasible - change bounds and resolve
2459    */
2460    if (newUpper[i]<1.0e100) {
2461      if(newLower[i]<1.0e100) {
2462        // feasible - no action
2463      } else {
2464        // up feasible, down infeasible
2465        returnCode=1;
2466        if (stopOnFirstInfeasible)
2467          break;
2468      }
2469    } else {
2470      if(newLower[i]<1.0e100) {
2471        // down feasible, up infeasible
2472        returnCode=1;
2473        if (stopOnFirstInfeasible)
2474          break;
2475      } else {
2476        // neither side feasible
2477        returnCode=-1;
2478        break;
2479      }
2480    }
2481  }
2482  delete [] rowSolution;
2483  delete [] columnSolution;
2484  delete [] saveStatus;
2485#endif
2486  objectiveValue_ = saveObjectiveValue;
2487  return returnCode;
2488}
2489// treat no pivot as finished (unless interesting)
2490int ClpSimplexDual::fastDual(bool alwaysFinish)
2491{
2492  algorithm_ = -1;
2493  dualTolerance_=dblParam_[ClpDualTolerance];
2494  primalTolerance_=dblParam_[ClpPrimalTolerance];
2495
2496  // save dual bound
2497  double saveDualBound = dualBound_;
2498
2499  double objectiveChange;
2500  // for dual we will change bounds using dualBound_
2501  // for this we need clean basis so it is after factorize
2502  gutsOfSolution(rowActivityWork_,columnActivityWork_);
2503  numberFake_ =0; // Number of variables at fake bounds
2504  changeBounds(true,NULL,objectiveChange);
2505
2506  problemStatus_ = -1;
2507  numberIterations_=0;
2508
2509  int lastCleaned=0; // last time objective or bounds cleaned up
2510
2511  // number of times we have declared optimality
2512  numberTimesOptimal_=0;
2513
2514  // Progress indicator
2515  ClpSimplexProgress progress(this);
2516
2517  // This says whether to restore things etc
2518  int factorType=0;
2519  /*
2520    Status of problem:
2521    0 - optimal
2522    1 - infeasible
2523    2 - unbounded
2524    -1 - iterating
2525    -2 - factorization wanted
2526    -3 - redo checking without factorization
2527    -4 - looks infeasible
2528
2529    BUT also from whileIterating return code is:
2530
2531   -1 iterations etc
2532   -2 inaccuracy
2533   -3 slight inaccuracy (and done iterations)
2534   +0 looks optimal (might be unbounded - but we will investigate)
2535   +1 looks infeasible
2536   +3 max iterations
2537
2538  */
2539
2540  int returnCode = 0;
2541
2542  while (problemStatus_<0) {
2543    int iRow,iColumn;
2544    // clear
2545    for (iRow=0;iRow<4;iRow++) {
2546      rowArray_[iRow]->clear();
2547    }   
2548   
2549    for (iColumn=0;iColumn<2;iColumn++) {
2550      columnArray_[iColumn]->clear();
2551    }   
2552
2553    // give matrix (and model costs and bounds a chance to be
2554    // refreshed (normally null)
2555    matrix_->refresh(this);
2556    // may factorize, checks if problem finished
2557    // should be able to speed this up on first time
2558    statusOfProblemInDual(lastCleaned,factorType,progress);
2559
2560    // Say good factorization
2561    factorType=1;
2562
2563    // Do iterations
2564    if (problemStatus_<0) {
2565      returnCode = whileIterating();
2566      if (!alwaysFinish&&returnCode<1) {
2567        double limit = 0.0;
2568        getDblParam(ClpDualObjectiveLimit, limit);
2569        if(objectiveValue()*optimizationDirection_<limit|| 
2570           numberAtFakeBound()) {
2571          // can't say anything interesting - might as well return
2572#ifdef CLP_DEBUG
2573          printf("returning from fastDual after %d iterations with code %d\n",
2574                 numberIterations_,returnCode);
2575#endif
2576          returnCode=1;
2577          break;
2578        }
2579      }
2580      returnCode=0;
2581    }
2582  }
2583
2584  assert(!numberFake_||returnCode||problemStatus_); // all bounds should be okay
2585  dualBound_ = saveDualBound;
2586  return returnCode;
2587}
2588/* Checks number of variables at fake bounds.  This is used by fastDual
2589   so can exit gracefully before end */
2590int 
2591ClpSimplexDual::numberAtFakeBound()
2592{
2593  int iSequence;
2594  int numberFake=0;
2595 
2596  for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
2597    FakeBound bound = getFakeBound(iSequence);
2598    switch(getStatus(iSequence)) {
2599
2600    case basic:
2601      break;
2602    case isFree:
2603    case superBasic:
2604      break;
2605    case atUpperBound:
2606      if (bound==upperFake||bound==bothFake)
2607        numberFake++;
2608      break;
2609    case atLowerBound:
2610      if (bound==lowerFake||bound==bothFake)
2611        numberFake++;
2612      break;
2613    }
2614  }
2615  numberFake_ = numberFake;
2616  return numberFake;
2617}
Note: See TracBrowser for help on using the repository browser.