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

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

Improve speed

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 64.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 "OsiPackedMatrix.hpp"
100#include "OsiIndexedVector.hpp"
101#include "OsiWarmStartBasis.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// This returns a non const array filled with input from scalar
110// or actual array
111template <class T> inline T*
112copyOfArray( const T * array, const int size, T value)
113{
114  T * arrayNew = new T[size];
115  if (array) 
116    CoinDisjointCopyN(array,size,arrayNew);
117  else
118    CoinFillN ( arrayNew, size,value);
119  return arrayNew;
120}
121
122// This returns a non const array filled with actual array (or NULL)
123template <class T> inline T*
124copyOfArray( const T * array, const int size)
125{
126  if (array) {
127    T * arrayNew = new T[size];
128    CoinDisjointCopyN(array,size,arrayNew);
129    return arrayNew;
130  } else {
131    return NULL;
132  }
133}
134// dual
135int ClpSimplexDual::dual ( )
136{
137
138  /* *** Method
139     This is a vanilla version of dual simplex.
140
141     It tries to be a single phase approach with a weight of 1.0 being
142     given to getting optimal and a weight of dualBound_ being
143     given to getting dual feasible.  In this version I have used the
144     idea that this weight can be thought of as a fake bound.  If the
145     distance between the lower and upper bounds on a variable is less
146     than the feasibility weight then we are always better off flipping
147     to other bound to make dual feasible.  If the distance is greater
148     then we make up a fake bound dualBound_ away from one bound.
149     If we end up optimal or primal infeasible, we check to see if
150     bounds okay.  If so we have finished, if not we increase dualBound_
151     and continue (after checking if unbounded). I am undecided about
152     free variables - there is coding but I am not sure about it.  At
153     present I put them in basis anyway.
154
155     The code is designed to take advantage of sparsity so arrays are
156     seldom zeroed out from scratch or gone over in their entirety.
157     The only exception is a full scan to find outgoing variable.  This
158     will be changed to keep an updated list of infeasibilities (or squares
159     if steepest edge).  Also on easy problems we don't need full scan - just
160     pick first reasonable.
161
162     One problem is how to tackle degeneracy and accuracy.  At present
163     I am using the modification of costs which I put in OSL and which was
164     extended by Gill et al.  I am still not sure of the exact details.
165
166     The flow of dual is three while loops as follows:
167
168     while (not finished) {
169
170       while (not clean solution) {
171
172          Factorize and/or clean up solution by flipping variables so
173          dual feasible.  If looks finished check fake dual bounds.
174          Repeat until status is iterating (-1) or finished (0,1,2)
175
176       }
177
178       while (status==-1) {
179
180         Iterate until no pivot in or out or time to re-factorize.
181
182         Flow is:
183
184         choose pivot row (outgoing variable).  if none then
185         we are primal feasible so looks as if done but we need to
186         break and check bounds etc.
187
188         Get pivot row in tableau
189
190         Choose incoming column.  If we don't find one then we look
191         primal infeasible so break and check bounds etc.  (Also the
192         pivot tolerance is larger after any iterations so that may be
193         reason)
194
195         If we do find incoming column, we may have to adjust costs to
196         keep going forwards (anti-degeneracy).  Check pivot will be stable
197         and if unstable throw away iteration (we will need to implement
198         flagging of basic variables sometime) and break to re-factorize.
199         If minor error re-factorize after iteration.
200
201         Update everything (this may involve flipping variables to stay
202         dual feasible.
203
204       }
205
206     }
207
208     At present we never check we are going forwards.  I overdid that in
209     OSL so will try and make a last resort.
210
211     Needs partial scan pivot out option.
212     Needs dantzig, uninitialized and full steepest edge options (can still
213     use partial scan)
214
215     May need other anti-degeneracy measures, especially if we try and use
216     loose tolerances as a way to solve in fewer iterations.
217
218     I like idea of dynamic scaling.  This gives opportunity to decouple
219     different implications of scaling for accuracy, iteration count and
220     feasibility tolerance.
221
222  */
223
224
225  // sanity check
226  assert (numberRows_==matrix_->getNumRows());
227  assert (numberColumns_==matrix_->getNumCols());
228  // for moment all arrays must exist
229  assert(columnLower_);
230  assert(columnUpper_);
231  assert(rowLower_);
232  assert(rowUpper_);
233
234#ifdef CLP_DEBUG
235  int debugIteration=-1;
236#endif
237
238  algorithm_ = -1;
239  dualTolerance_=dblParam_[OsiDualTolerance];
240  primalTolerance_=dblParam_[OsiPrimalTolerance];
241
242  // put in standard form (and make row copy)
243  // create modifiable copies of model rim and do optional scaling
244  createRim(7+8+16,true);
245
246  // save dual bound
247  double saveDualBound_ = dualBound_;
248
249  int iRow,iColumn;
250  // Do initial factorization
251  // and set certain stuff
252  // We can either set increasing rows so ...IsBasic gives pivot row
253  // or we can just increment iBasic one by one
254  // for now let ...iBasic give pivot row
255  factorization_->increasingRows(2);
256  // row activities have negative sign
257  factorization_->slackValue(-1.0);
258  factorization_->zeroTolerance(1.0e-13);
259  // save if sparse factorization wanted
260  int saveSparse = factorization_->sparseThreshold();
261
262  if (internalFactorize(0))
263    return 1; // some error
264
265  // If user asked for perturbation - do it
266  int savePerturbation = perturbation_;
267
268  if (perturbation_<100) 
269    perturb();
270
271  double objectiveChange;
272  // for dual we will change bounds using dualBound_
273  // for this we need clean basis so it is after factorize
274  gutsOfSolution(rowActivityWork_,columnActivityWork_);
275  changeBounds(true,NULL,objectiveChange);
276
277  problemStatus_ = -1;
278  numberIterations_=0;
279
280  int lastCleaned=0; // last time objective or bounds cleaned up
281
282  // number of times we have declared optimality
283  numberTimesOptimal_=0;
284
285  // This says whether to restore things etc
286  int factorType=0;
287  /*
288    Status of problem:
289    0 - optimal
290    1 - infeasible
291    2 - unbounded
292    -1 - iterating
293    -2 - factorization wanted
294    -3 - redo checking without factorization
295    -4 - looks infeasible
296  */
297  while (problemStatus_<0) {
298    // clear
299    for (iRow=0;iRow<4;iRow++) {
300      rowArray_[iRow]->clear();
301    }   
302   
303    for (iColumn=0;iColumn<2;iColumn++) {
304      columnArray_[iColumn]->clear();
305    }   
306
307    // give matrix (and model costs and bounds a chance to be
308    // refreshed (normally null)
309    matrix_->refresh(this);
310    // If getting nowhere - why not give it a kick
311#if 0
312    // does not seem to work too well - do some more work
313    if (perturbation_<101&&numberIterations_>2*(numberRows_+numberColumns_))
314      perturb();
315#endif
316    // may factorize, checks if problem finished
317    statusOfProblemInDual(lastCleaned,factorType);
318
319    // Say good factorization
320    factorType=1;
321    if (saveSparse) {
322      // use default at present
323      factorization_->sparseThreshold(0);
324      factorization_->goSparse();
325    }
326
327    // status stays at -1 while iterating, >=0 finished, -2 to invert
328    // status -3 to go to top without an invert
329    while (problemStatus_==-1) {
330#ifdef CLP_DEBUG
331      {
332        int i;
333        for (i=0;i<4;i++) {
334          rowArray_[i]->checkClear();
335        }   
336        for (i=0;i<2;i++) {
337          columnArray_[i]->checkClear();
338        }   
339      }     
340#endif
341#if CLP_DEBUG>2
342      // very expensive
343      if (numberIterations_>0&&numberIterations_<-801) {
344        handler_->setLogLevel(63);
345        double saveValue = objectiveValue_;
346        double * saveRow1 = new double[numberRows_];
347        double * saveRow2 = new double[numberRows_];
348        memcpy(saveRow1,rowReducedCost_,numberRows_*sizeof(double));
349        memcpy(saveRow2,rowActivityWork_,numberRows_*sizeof(double));
350        double * saveColumn1 = new double[numberColumns_];
351        double * saveColumn2 = new double[numberColumns_];
352        memcpy(saveColumn1,reducedCostWork_,numberColumns_*sizeof(double));
353        memcpy(saveColumn2,columnActivityWork_,numberColumns_*sizeof(double));
354        gutsOfSolution(rowActivityWork_,columnActivityWork_);
355        printf("xxx %d old obj %g, recomputed %g, sum dual inf %g\n",
356               numberIterations_,
357               saveValue,objectiveValue_,sumDualInfeasibilities_);
358        if (saveValue>objectiveValue_+1.0e-2)
359          printf("**bad**\n");
360        memcpy(rowReducedCost_,saveRow1,numberRows_*sizeof(double));
361        memcpy(rowActivityWork_,saveRow2,numberRows_*sizeof(double));
362        memcpy(reducedCostWork_,saveColumn1,numberColumns_*sizeof(double));
363        memcpy(columnActivityWork_,saveColumn2,numberColumns_*sizeof(double));
364        delete [] saveRow1;
365        delete [] saveRow2;
366        delete [] saveColumn1;
367        delete [] saveColumn2;
368        objectiveValue_=saveValue;
369      }
370#endif
371#ifdef CLP_DEBUG
372      {
373        int iSequence, number=numberRows_+numberColumns_;
374        for (iSequence=0;iSequence<number;iSequence++) {
375          double lowerValue=lower_[iSequence];
376          double upperValue=upper_[iSequence];
377          double value=solution_[iSequence];
378          if(getStatus(iSequence)!=ClpSimplex::basic) {
379            assert(lowerValue>-1.0e20);
380            assert(upperValue<1.0e20);
381          }
382          switch(getStatus(iSequence)) {
383           
384          case ClpSimplex::basic:
385            break;
386          case ClpSimplex::isFree:
387          case ClpSimplex::superBasic:
388            break;
389          case ClpSimplex::atUpperBound:
390            assert (fabs(value-upperValue)<=primalTolerance_) ;
391            break;
392          case ClpSimplex::atLowerBound:
393            assert (fabs(value-lowerValue)<=primalTolerance_) ;
394              break;
395          }
396        }
397      }
398      if(numberIterations_==debugIteration) {
399        printf("dodgy iteration coming up\n");
400      }
401#endif
402      // choose row to go out
403      dualRow();
404      if (pivotRow_>=0) {
405        // we found a pivot row
406        handler_->message(CLP_SIMPLEX_PIVOTROW,messages_)
407          <<pivotRow_
408          <<OsiMessageEol;
409        // check accuracy of weights
410        dualRowPivot_->checkAccuracy();
411        // get sign for finding row of tableau
412        rowArray_[0]->insert(pivotRow_,directionOut_);
413        factorization_->updateColumnTranspose(rowArray_[1],rowArray_[0]);
414        // put row of tableau in rowArray[0] and columnArray[0]
415        matrix_->transposeTimes(this,-1.0,
416                                rowArray_[0],columnArray_[1],columnArray_[0]);
417        // rowArray has pi equivalent
418        // do ratio test
419        dualColumn(rowArray_[0],columnArray_[0],columnArray_[1],
420                   rowArray_[3]);
421        if (sequenceIn_>=0) {
422          // normal iteration
423          // update the incoming column
424          unpack(rowArray_[1]);
425          factorization_->updateColumn(rowArray_[2],rowArray_[1],true);
426          // and update dual weights (can do in parallel - with extra array)
427          dualRowPivot_->updateWeights(rowArray_[0],rowArray_[2],
428                                       rowArray_[1]);
429          // see if update stable
430          double btranAlpha = -alpha_*directionOut_; // for check
431          alpha_=(*rowArray_[1])[pivotRow_];
432#ifdef CLP_DEBUG
433          if ((handler_->logLevel()&32))
434            printf("btran alpha %g, ftran alpha %g\n",btranAlpha,alpha_);
435#endif
436          if (fabs(btranAlpha)<1.0e-12||fabs(alpha_)<1.0e-12||
437              fabs(btranAlpha-alpha_)>1.0e-7*(1.0+fabs(alpha_))) {
438            handler_->message(CLP_DUAL_CHECK,messages_)
439              <<btranAlpha
440              <<alpha_
441              <<OsiMessageEol;
442            dualRowPivot_->unrollWeights();
443            if (factorization_->pivots()) {
444              problemStatus_=-2; // factorize now
445              rowArray_[0]->clear();
446              rowArray_[1]->clear();
447              columnArray_[0]->clear();
448              break;
449            } else {
450              // take on more relaxed criterion
451              if (fabs(btranAlpha)<1.0e-12||fabs(alpha_)<1.0e-12||
452                  fabs(btranAlpha-alpha_)>1.0e-4*(1.0+fabs(alpha_))) {
453                // need to reject something
454                char x = isColumn(sequenceOut_) ? 'C' :'R';
455                handler_->message(CLP_SIMPLEX_FLAG,messages_)
456                  <<x<<sequenceWithin(sequenceOut_)
457                  <<OsiMessageEol;
458                setFlagged(sequenceOut_);
459                lastBadIteration_ = numberIterations_; // say be more cautious
460                rowArray_[0]->clear();
461                rowArray_[1]->clear();
462                columnArray_[0]->clear();
463                continue;
464              }
465            }
466          }
467          // update duals BEFORE replaceColumn so can do updateColumn
468          objectiveChange=0.0;
469          // do duals first as variables may flip bounds
470          // rowArray_[0] and columnArray_[0] may have flips
471          // so use rowArray_[3] for work array from here on
472          int nswapped = 
473            updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[2],theta_,
474                        objectiveChange);
475          // which will change basic solution
476          if (nswapped) {
477#ifdef CLP_DEBUG
478          if ((handler_->logLevel()&16))
479            printf("old dualOut_ %g, v %g, l %g, u %g - new ",
480                   dualOut_,valueOut_,lowerOut_,upperOut_);
481            double oldOut=dualOut_;
482#endif
483            factorization_->updateColumn(rowArray_[3],rowArray_[2],false);
484            dualRowPivot_->updatePrimalSolution(rowArray_[2],
485                                                1.0,objectiveChange);
486
487            // recompute dualOut_
488            valueOut_ = solution_[sequenceOut_];
489            if (directionOut_<0) {
490              dualOut_ = valueOut_ - upperOut_;
491            } else {
492              dualOut_ = lowerOut_ - valueOut_;
493            }
494#ifdef CLP_DEBUG
495            if ((handler_->logLevel()&16))
496              printf("%g\n",dualOut_);
497            assert(dualOut_<=oldOut);
498#endif
499            if(dualOut_<0.0&&factorization_->pivots()) {
500              // going backwards - factorize
501              dualRowPivot_->unrollWeights();
502              problemStatus_=-2; // factorize now
503              break;
504            }
505          }
506          // amount primal will move
507          double movement = -dualOut_*directionOut_/alpha_;
508          // so objective should increase by fabs(dj)*movement
509          // but we already have objective change - so check will be good
510          if (objectiveChange+fabs(movement*dualIn_)<-1.0e-5) {
511#ifdef CLP_DEBUG
512            if (handler_->logLevel()&32)
513              printf("movement %g, swap change %g, rest %g  * %g\n",
514                     objectiveChange+fabs(movement*dualIn_),
515                     objectiveChange,movement,dualIn_);
516#endif
517            if(factorization_->pivots()>5) {
518              // going backwards - factorize
519              dualRowPivot_->unrollWeights();
520              problemStatus_=-2; // factorize now
521              break;
522            }
523          }
524          // if stable replace in basis
525          int updateStatus = factorization_->replaceColumn(rowArray_[2],
526                                                           pivotRow_,
527                                                           alpha_);
528          if (updateStatus==1) {
529            // slight error
530            if (factorization_->pivots()>5)
531              problemStatus_=-2; // factorize now
532          } else if (updateStatus==2) {
533            // major error
534            dualRowPivot_->unrollWeights();
535            // later we may need to unwind more e.g. fake bounds
536            if (factorization_->pivots()) {
537              problemStatus_=-2; // factorize now
538              break;
539            } else {
540              // need to reject something
541              char x = isColumn(sequenceOut_) ? 'C' :'R';
542              handler_->message(CLP_SIMPLEX_FLAG,messages_)
543                <<x<<sequenceWithin(sequenceOut_)
544                <<OsiMessageEol;
545              setFlagged(sequenceOut_);
546              lastBadIteration_ = numberIterations_; // say be more cautious
547              rowArray_[0]->clear();
548              rowArray_[1]->clear();
549              columnArray_[0]->clear();
550              continue;
551            }
552          } else if (updateStatus==3) {
553            // out of memory
554            // increase space if not many iterations
555            if (factorization_->pivots()<
556                0.5*factorization_->maximumPivots()&&
557                factorization_->pivots()<200)
558              factorization_->areaFactor(
559                                         factorization_->areaFactor() * 1.1);
560            problemStatus_=-2; // factorize now
561          } 
562          // update primal solution
563          if (theta_<0.0) {
564#ifdef CLP_DEBUG
565            if (handler_->logLevel()&32)
566              printf("negative theta %g\n",theta_);
567#endif
568            theta_=0.0;
569          }
570          // do actual flips
571          flipBounds(rowArray_[0],columnArray_[0],theta_);
572          dualRowPivot_->updatePrimalSolution(rowArray_[1],
573                                              movement,
574                                              objectiveChange);
575#ifdef CLP_DEBUG
576          double oldobj=objectiveValue_;
577#endif
578          int whatNext=housekeeping(objectiveChange);
579          // and set bounds correctly
580          originalBound(sequenceIn_); 
581          changeBound(sequenceOut_);
582#ifdef CLP_DEBUG
583          if (objectiveValue_<oldobj-1.0e-5&&(handler_->logLevel()&16))
584            printf("obj backwards %g %g\n",objectiveValue_,oldobj);
585#endif
586          if (whatNext==1) {
587            problemStatus_ =-2; // refactorize
588          } else if (whatNext==2) {
589            // maximum iterations or equivalent
590            problemStatus_= 3;
591            break;
592          }
593        } else {
594          // no incoming column is valid
595#ifdef CLP_DEBUG
596          if (handler_->logLevel()&32)
597            printf("** no column pivot\n");
598#endif
599          if (factorization_->pivots()<5) {
600            problemStatus_=-4; //say looks infeasible
601            // create ray anyway
602            delete [] ray_;
603            ray_ = new double [ numberRows_];
604            CoinDisjointCopyN(rowArray_[0]->denseVector(),numberRows_,ray_);
605          }
606          rowArray_[0]->clear();
607          columnArray_[0]->clear();
608          break;
609        }
610      } else {
611        // no pivot row
612#ifdef CLP_DEBUG
613        if (handler_->logLevel()&32)
614          printf("** no row pivot\n");
615#endif
616        if (!factorization_->pivots()) {
617          // may have crept through - so may be optimal
618          //problemStatus_=-5; //say looks unbounded
619          problemStatus_=0;
620          // check any flagged variables
621          int iRow;
622          for (iRow=0;iRow<numberRows_;iRow++) {
623            int iPivot=pivotVariable_[iRow];
624            if (flagged(iPivot))
625              break;
626          }
627          if (iRow<numberRows_) {
628#ifdef CLP_DEBUG
629            std::cerr<<"Flagged variables at end - infeasible?"<<std::endl;
630#endif
631            problemStatus_=-4; //say looks infeasible
632            // create ray anyway
633            delete [] ray_;
634            ray_ = new double [ numberRows_];
635            CoinDisjointCopyN(rowArray_[0]->denseVector(),numberRows_,ray_);
636          }
637        }
638        break;
639      }
640    }
641  }
642
643  // at present we are leaving factorization around
644  // maybe we should empty it
645  deleteRim();
646  handler_->message(CLP_SIMPLEX_FINISHED+problemStatus_,messages_)
647    <<objectiveValue()
648    <<OsiMessageEol;
649  // Restore any saved stuff
650  perturbation_ = savePerturbation;
651  factorization_->sparseThreshold(saveSparse);
652  dualBound_ = saveDualBound_;
653  return problemStatus_;
654}
655/* The duals are updated by the given arrays.
656    Returns number of infeasibilities.
657    rowArray and columnarray will have flipped
658    The output vector has movement (row length array) */
659int
660ClpSimplexDual::updateDualsInDual(OsiIndexedVector * rowArray,
661                                OsiIndexedVector * columnArray,
662                                OsiIndexedVector * outputArray,
663                                double theta,
664                                double & objectiveChange)
665{
666
667  outputArray->clear();
668 
669  double * work;
670  int number;
671  int * which;
672 
673  int numberInfeasibilities=0;
674  int numberRowInfeasibilities=0;
675 
676  // see whether we will be doing full recompute
677  bool fullRecompute= (rowArray->getNumElements()==numberRows_&&
678                       columnArray->getNumElements()==numberColumns_);
679  int numberAtFake=0;
680
681  // use a tighter tolerance except for all being okay
682  double tolerance = dualTolerance_;
683
684  double changeObj=0.0;
685
686  int iSection;
687 
688  for (iSection=0;iSection<2;iSection++) {
689    int i;
690    double * solution = solutionRegion(iSection);
691    double * reducedCost = djRegion(iSection);
692    double * lower = lowerRegion(iSection);
693    double * upper = upperRegion(iSection);
694    double * cost = costRegion(iSection);
695    int addSequence;
696    if (!iSection) {
697      addSequence = numberColumns_;
698      work = rowArray->denseVector();
699      number = rowArray->getNumElements();
700      which = rowArray->getIndices();
701    } else {
702      // set number of infeasibilities in row array
703      addSequence=0;
704      numberRowInfeasibilities=numberInfeasibilities;
705      rowArray->setNumElements(numberInfeasibilities);
706      numberInfeasibilities=0;
707      work = columnArray->denseVector();
708      number = columnArray->getNumElements();
709      which = columnArray->getIndices();
710    }
711   
712    for (i=0;i<number;i++) {
713      int iSequence = which[i];
714      double alphaI = work[iSequence];
715      double value = reducedCost[iSequence]-theta*alphaI;
716      work[iSequence]=0.0;
717      reducedCost[iSequence]=value;
718     
719      if (!fixed(iSequence+addSequence)) {
720        double movement=0.0;
721        FakeBound bound = getFakeBound(iSequence+addSequence);
722        Status status = getStatus(iSequence+addSequence);
723       
724        switch(status) {
725         
726        case ClpSimplex::basic:
727        case ClpSimplex::superBasic:
728          break;
729        case ClpSimplex::isFree:
730          if (fabs(value)>tolerance) { 
731#ifdef CLP_DEBUG
732            if (handler_->logLevel()&32)
733              printf("%d %d, free has dj of %g, alpha %g\n",
734                     iSection,iSequence,value,alphaI);
735#endif
736          }
737          break;
738        case ClpSimplex::atUpperBound:
739          if (value>tolerance) {
740            // to lower bound (if swap)
741            // put back alpha
742            which[numberInfeasibilities++]=iSequence;
743            work[iSequence]=alphaI;
744            movement = lower[iSequence]-upper[iSequence];
745#ifdef CLP_DEBUG
746            if ((handler_->logLevel()&32))
747              printf("%d %d, new dj %g, alpha %g, movement %g\n",
748                   iSection,iSequence,value,alphaI,movement);
749#endif
750            changeObj += movement*cost[iSequence];
751            if (bound==ClpSimplexDual::bothFake||
752                bound==ClpSimplexDual::lowerFake) 
753              numberAtFake++;
754          } else if (fullRecompute) {
755            // at correct bound
756            if (bound==ClpSimplexDual::bothFake||
757                bound==ClpSimplexDual::upperFake) {
758              // but flip if dj would allow
759              if (bound==ClpSimplexDual::upperFake&&
760                  value>=-tolerance) {
761                movement = lower[iSequence]-upper[iSequence];
762                setStatus(iSequence+addSequence,ClpSimplex::atLowerBound);
763                solution[iSequence] = lower[iSequence];
764                changeObj += movement*cost[iSequence];
765              } else {
766                numberAtFake++;
767              }
768            }
769          }
770          break;
771        case ClpSimplex::atLowerBound:
772          if (value<-tolerance) {
773            // to upper bound
774            // put back alpha
775            which[numberInfeasibilities++]=iSequence;
776            work[iSequence]=alphaI;
777            movement = upper[iSequence] - lower[iSequence];
778#ifdef CLP_DEBUG
779            if ((handler_->logLevel()&32))
780              printf("%d %d, new dj %g, alpha %g, movement %g\n",
781                   iSection,iSequence,value,alphaI,movement);
782#endif
783            changeObj += movement*cost[iSequence];
784            if (bound==ClpSimplexDual::bothFake||
785                bound==ClpSimplexDual::upperFake) 
786              numberAtFake++;
787          } else if (fullRecompute) {
788            // at correct bound
789            if (bound==ClpSimplexDual::bothFake||
790                bound==ClpSimplexDual::lowerFake) {
791              // but flip if dj would allow
792              if (bound==ClpSimplexDual::lowerFake&&
793                  value<=tolerance) {
794                movement = upper[iSequence] - lower[iSequence];
795                setStatus(iSequence+addSequence,ClpSimplex::atUpperBound);
796                solution[iSequence] = upper[iSequence];
797                changeObj += movement*cost[iSequence];
798              } else {
799                numberAtFake++;
800              }
801            }
802          }
803          break;
804        }
805        if (!fullRecompute) {
806          if (movement) {
807            if (!iSection) {
808              // row (sign ?)
809              outputArray->quickAdd(iSequence,-movement);
810            } else {
811              matrix_->add(this,outputArray,iSequence,movement);
812            }
813          }
814        }
815      }
816    }
817  }
818#ifdef CLP_DEBUG
819  if (fullRecompute&&numberAtFake&&(handler_->logLevel()&16)!=0) 
820    printf("%d fake after full update\n",numberAtFake);
821#endif
822  outputArray->stopQuickAdd();
823  // set number of infeasibilities
824  columnArray->setNumElements(numberInfeasibilities);
825  numberInfeasibilities += numberRowInfeasibilities;
826  if (fullRecompute) {
827    // do actual flips
828    flipBounds(rowArray,columnArray,theta);
829  }
830  objectiveChange += changeObj;
831  return numberInfeasibilities;
832}
833/*
834   Chooses dual pivot row
835   Would be faster with separate region to scan
836   and will have this (with square of infeasibility) when steepest
837   For easy problems we can just choose one of the first rows we look at
838*/
839void
840ClpSimplexDual::dualRow()
841{
842  // get pivot row using whichever method it is
843  pivotRow_=dualRowPivot_->pivotRow();
844  if (pivotRow_>=0) {
845    int iPivot=pivotVariable_[pivotRow_];
846    sequenceOut_ = iPivot;
847    if (iPivot>=numberColumns_) {
848      // slack
849      iPivot-=numberColumns_;
850      valueOut_=rowActivityWork_[iPivot];
851      lowerOut_=rowLowerWork_[iPivot];
852      upperOut_=rowUpperWork_[iPivot];
853    } else {
854      // column
855      valueOut_=columnActivityWork_[iPivot];
856      lowerOut_=columnLowerWork_[iPivot];
857      upperOut_=columnUpperWork_[iPivot];
858    }
859    // if we have problems we could try other way and hope we get a
860    // zero pivot?
861    if (valueOut_>upperOut_) {
862      directionOut_ = -1;
863      dualOut_ = valueOut_ - upperOut_;
864    } else {
865      directionOut_ = 1;
866      dualOut_ = lowerOut_ - valueOut_;
867    }
868#ifdef CLP_DEBUG
869    assert(dualOut_>=0.0);
870#endif
871  }
872  return ;
873}
874// Checks if any fake bounds active - if so returns number and modifies
875// dualBound_ and everything.
876// Free variables will be left as free
877// Returns number of bounds changed if >=0
878// Returns -1 if not initialize and no effect
879// Fills in changeVector which can be used to see if unbounded
880// and cost of change vector
881int
882ClpSimplexDual::changeBounds(bool initialize,
883                                 OsiIndexedVector * outputArray,
884                                 double & changeCost)
885{
886  if (!initialize) {
887    int numberInfeasibilities;
888    double newBound;
889    newBound = 5.0*dualBound_;
890    numberInfeasibilities=0;
891    changeCost=0.0;
892    // put back original bounds and then check
893    createRim(3);
894    int iSequence;
895    // bounds will get bigger - just look at ones at bounds
896    for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
897      double lowerValue=lower_[iSequence];
898      double upperValue=upper_[iSequence];
899      double value=solution_[iSequence];
900      setFakeBound(iSequence,ClpSimplexDual::noFake);
901      switch(getStatus(iSequence)) {
902       
903      case ClpSimplex::basic:
904        break;
905      case ClpSimplex::isFree:
906      case ClpSimplex::superBasic:
907        break;
908      case ClpSimplex::atUpperBound:
909        if (fabs(value-upperValue)>primalTolerance_) 
910          numberInfeasibilities++;
911        break;
912      case ClpSimplex::atLowerBound:
913        if (fabs(value-lowerValue)>primalTolerance_) 
914          numberInfeasibilities++;
915        break;
916      }
917    }
918    if (numberInfeasibilities) {
919      int iSequence;
920      for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
921        double lowerValue=lower_[iSequence];
922        double upperValue=upper_[iSequence];
923        double newLowerValue;
924        double newUpperValue;
925        Status status = getStatus(iSequence);
926        if (status==ClpSimplex::atUpperBound||
927            status==ClpSimplex::atLowerBound) {
928          double value = solution_[iSequence];
929          if (value-lowerValue<=upperValue-value) {
930            newLowerValue = max(lowerValue,value-0.666667*newBound);
931            newUpperValue = min(upperValue,newLowerValue+newBound);
932          } else {
933            newUpperValue = min(upperValue,value+0.666667*newBound);
934            newLowerValue = max(lowerValue,newUpperValue-newBound);
935          }
936          lower_[iSequence]=newLowerValue;
937          upper_[iSequence]=newUpperValue;
938          if (newLowerValue > lowerValue) {
939            if (newUpperValue < upperValue) 
940              setFakeBound(iSequence,ClpSimplexDual::bothFake);
941            else
942              setFakeBound(iSequence,ClpSimplexDual::lowerFake);
943          } else {
944            if (newUpperValue < upperValue) 
945              setFakeBound(iSequence,ClpSimplexDual::upperFake);
946          }
947          if (status==ClpSimplex::atUpperBound)
948            solution_[iSequence] = newUpperValue;
949          else
950            solution_[iSequence] = newLowerValue;
951          double movement = solution_[iSequence] - value;
952          if (movement&&outputArray) {
953            if (iSequence>=numberColumns_) {
954              outputArray->quickAdd(iSequence,-movement);
955              changeCost += movement*cost_[iSequence];
956            } else {
957              matrix_->add(this,outputArray,iSequence,movement);
958              changeCost += movement*cost_[iSequence];
959            }
960          }
961        }
962      }
963      dualBound_ = newBound;
964      if (outputArray)
965        outputArray->stopQuickAdd();
966    } else {
967      numberInfeasibilities=-1;
968    }
969    return numberInfeasibilities;
970  } else {
971    int iSequence;
972     
973    for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
974      Status status = getStatus(iSequence);
975      if (status==ClpSimplex::atUpperBound||
976          status==ClpSimplex::atLowerBound) {
977        double lowerValue=lower_[iSequence];
978        double upperValue=upper_[iSequence];
979        double value = solution_[iSequence];
980        if (lowerValue>-largeValue_||upperValue<largeValue_) {
981          if (lowerValue-value>-0.5*dualBound_||
982              upperValue-value<0.5*dualBound_) {
983            if (fabs(lowerValue-value)<=fabs(upperValue-value)) {
984              if (upperValue > lowerValue + dualBound_) {
985                upper_[iSequence]=lowerValue+dualBound_;
986                setFakeBound(iSequence,ClpSimplexDual::upperFake);
987              }
988            } else {
989              if (lowerValue < upperValue - dualBound_) {
990                lower_[iSequence]=upperValue-dualBound_;
991                setFakeBound(iSequence,ClpSimplexDual::lowerFake);
992              }
993            }
994          } else {
995            lower_[iSequence]=-0.5*dualBound_;
996            upper_[iSequence]= 0.5*dualBound_;
997            setFakeBound(iSequence,ClpSimplexDual::bothFake);
998          }
999        }
1000      }
1001    }
1002    return 1;
1003  }
1004}
1005/*
1006   Row array has row part of pivot row (as duals so sign may be switched)
1007   Column array has column part.
1008   This chooses pivot column.
1009   Spare array will be needed when we start getting clever.
1010   We will check for basic so spare array will never overflow.
1011   If necessary will modify costs
1012*/
1013void
1014ClpSimplexDual::dualColumn(OsiIndexedVector * rowArray,
1015                       OsiIndexedVector * columnArray,
1016                       OsiIndexedVector * spareArray,
1017                       OsiIndexedVector * spareArray2)
1018{
1019  double * work;
1020  int number;
1021  int * which;
1022  double * reducedCost;
1023
1024  int iSection;
1025
1026  sequenceIn_=-1;
1027  int numberPossiblySwapped=0;
1028  int numberRemaining=0;
1029 
1030  double totalThru=0.0; // for when variables flip
1031  double acceptablePivot=1.0e-7;
1032  if (factorization_->pivots())
1033    acceptablePivot=1.0e-5; // if we have iterated be more strict
1034  double bestEverPivot=acceptablePivot;
1035  int lastSequence = -1;
1036  double lastPivot=0.0;
1037  double upperTheta;
1038  double newTolerance = dualTolerance_;
1039  // will we need to increase tolerance
1040  bool thisIncrease=false;
1041  // If we think we need to modify costs (not if something from broad sweep)
1042  bool modifyCosts=false;
1043  // Increase in objective due to swapping bounds (may be negative)
1044  double increaseInObjective=0.0;
1045
1046  // use spareArrays to put ones looked at in
1047  // we are going to flip flop between
1048  int iFlip = 0;
1049  // Possible list of pivots
1050  int interesting[2];
1051  // where possible swapped ones are
1052  int swapped[2];
1053  // for zeroing out arrays after
1054  int marker[2][2];
1055  // pivot elements
1056  double * array[2], * spare, * spare2;
1057  // indices
1058  int * indices[2], * index, * index2;
1059  spareArray->clear();
1060  spareArray2->clear();
1061  array[0] = spareArray->denseVector();
1062  indices[0] = spareArray->getIndices();
1063  spare = array[0];
1064  index = indices[0];
1065  array[1] = spareArray2->denseVector();
1066  indices[1] = spareArray2->getIndices();
1067  int i;
1068  double * lower;
1069  double * upper;
1070
1071  // initialize lists
1072  for (i=0;i<2;i++) {
1073    interesting[i]=0;
1074    swapped[i]=numberColumns_;
1075    marker[i][0]=0;
1076    marker[i][1]=numberColumns_;
1077  }
1078
1079  /*
1080    First we get a list of possible pivots.  We can also see if the
1081    problem looks infeasible or whether we want to pivot in free variable.
1082    This may make objective go backwards but can only happen a finite
1083    number of times and I do want free variables basic.
1084
1085    Then we flip back and forth.  At the start of each iteration
1086    interesting[iFlip] should have possible candidates and swapped[iFlip]
1087    will have pivots if we decide to take a previous pivot.
1088    At end of each iteration interesting[1-iFlip] should have
1089    candidates if we go through this theta and swapped[1-iFlip]
1090    pivots if we don't go through.
1091
1092    At first we increase theta and see what happens.  We start
1093    theta at a reasonable guess.  If in right area then we do bit by bit.
1094
1095   */
1096
1097  // do first pass to get possibles
1098  // We can also see if infeasible or pivoting on free
1099  double tentativeTheta = 1.0e22;
1100  upperTheta = 1.0e31;
1101  double freePivot = acceptablePivot;
1102  for (iSection=0;iSection<2;iSection++) {
1103
1104    int addSequence;
1105
1106    if (!iSection) {
1107      lower = rowLowerWork_;
1108      upper = rowUpperWork_;
1109      work = rowArray->denseVector();
1110      number = rowArray->getNumElements();
1111      which = rowArray->getIndices();
1112      reducedCost = rowReducedCost_;
1113      addSequence = numberColumns_;
1114    } else {
1115      lower = columnLowerWork_;
1116      upper = columnUpperWork_;
1117      work = columnArray->denseVector();
1118      number = columnArray->getNumElements();
1119      which = columnArray->getIndices();
1120      reducedCost = reducedCostWork_;
1121      addSequence = 0;
1122    }
1123   
1124    for (i=0;i<number;i++) {
1125      int iSequence = which[i];
1126      double alpha = work[iSequence];
1127      if (fixed(iSequence+addSequence)||!alpha) 
1128        continue; // skip fixed ones or (zeroed out)
1129      double oldValue = reducedCost[iSequence];
1130      double value = oldValue-tentativeTheta*alpha;
1131      int keep = 0;
1132       
1133      switch(getStatus(iSequence+addSequence)) {
1134         
1135      case ClpSimplex::basic:
1136        break;
1137      case ClpSimplex::isFree:
1138      case ClpSimplex::superBasic:
1139        if (oldValue>dualTolerance_) {
1140          if (value<-newTolerance) 
1141            keep = 2;
1142        } else if (oldValue<-dualTolerance_) {
1143          if (value>newTolerance) 
1144            keep = 2;
1145        } else {
1146          if (alpha>=acceptablePivot) 
1147            keep = 2;
1148          else if (-alpha>=acceptablePivot) 
1149            keep = 2;
1150        }
1151        break;
1152      case ClpSimplex::atUpperBound:
1153        assert (oldValue<=dualTolerance_*1.0001);
1154        if (value>newTolerance) {
1155          keep = 1;
1156          value = oldValue-upperTheta*alpha;
1157          if (value>newTolerance && -alpha>=acceptablePivot) 
1158            upperTheta = (oldValue-newTolerance)/alpha;
1159        }
1160        break;
1161      case ClpSimplex::atLowerBound:
1162        assert (oldValue>=-dualTolerance_*1.0001);
1163        if (value<-newTolerance) {
1164          keep = 1;
1165          value = oldValue-upperTheta*alpha;
1166          if (value<-newTolerance && alpha>=acceptablePivot) 
1167            upperTheta = (oldValue+newTolerance)/alpha;
1168        }
1169        break;
1170      }
1171      if (keep) {
1172        if (keep==2) {
1173          // free - choose largest
1174          if (fabs(alpha)>freePivot) {
1175            freePivot=fabs(alpha);
1176            sequenceIn_ = iSequence + addSequence;
1177            theta_=oldValue/alpha;
1178          }
1179        } else {
1180          // add to list
1181          spare[numberRemaining]=alpha;
1182          index[numberRemaining++]=iSequence+addSequence;
1183        }
1184      }
1185    }
1186  }
1187  interesting[0]=numberRemaining;
1188  marker[0][0] = numberRemaining;
1189
1190  if (!numberRemaining)
1191    return; // Looks infeasible
1192
1193  if (sequenceIn_>=0) {
1194    // free variable - always choose
1195  } else {
1196
1197    theta_=1.0e50;
1198    // now flip flop between spare arrays until reasonable theta
1199    tentativeTheta = max(10.0*upperTheta,1.0e-7);
1200   
1201    // loops increasing tentative theta until can't go through
1202   
1203    while (tentativeTheta < 1.0e22) {
1204      double thruThis = 0.0;
1205     
1206      double bestPivot=acceptablePivot;
1207      int bestSequence=-1;
1208     
1209      numberPossiblySwapped = numberColumns_;
1210      numberRemaining = 0;
1211
1212      upperTheta = 1.0e50;
1213
1214      spare = array[iFlip];
1215      index = indices[iFlip];
1216      spare2 = array[1-iFlip];
1217      index2 = indices[1-iFlip];
1218
1219      // try 3 different ways
1220      // 1 bias increase by ones with slightly wrong djs
1221      // 2 bias by all
1222      // 3 bias by all - tolerance (doesn't seem very good)
1223#define TRYBIAS 1
1224
1225
1226      double increaseInThis=0.0; //objective increase in this loop
1227     
1228      for (i=0;i<interesting[iFlip];i++) {
1229        int iSequence = index[i];
1230        double alpha = spare[i];
1231        double oldValue = dj_[iSequence];
1232        double value = oldValue-tentativeTheta*alpha;
1233
1234        if (alpha < 0.0) {
1235          //at upper bound
1236          if (value>newTolerance) {
1237            double range = upper_[iSequence] - lower_[iSequence];
1238            thruThis -= range*alpha;
1239#if TRYBIAS==1
1240            if (oldValue>0.0)
1241              increaseInThis -= oldValue*range;
1242#elif TRYBIAS==2
1243            increaseInThis -= oldValue*range;
1244#else
1245            increaseInThis -= (oldValue+dualTolerance_)*range;
1246#endif
1247            // goes on swapped list (also means candidates if too many)
1248            spare2[--numberPossiblySwapped]=alpha;
1249            index2[numberPossiblySwapped]=iSequence;
1250            if (fabs(alpha)>bestPivot) {
1251              bestPivot=fabs(alpha);
1252              bestSequence=numberPossiblySwapped;
1253            }
1254          } else {
1255            value = oldValue-upperTheta*alpha;
1256            if (value>newTolerance && -alpha>=acceptablePivot) 
1257              upperTheta = (oldValue-newTolerance)/alpha;
1258            spare2[numberRemaining]=alpha;
1259            index2[numberRemaining++]=iSequence;
1260          }
1261        } else {
1262          // at lower bound
1263          if (value<-newTolerance) {
1264            double range = upper_[iSequence] - lower_[iSequence];
1265            thruThis += range*alpha;
1266            //?? is this correct - and should we look at good ones
1267#if TRYBIAS==1
1268            if (oldValue<0.0)
1269              increaseInThis += oldValue*range;
1270#elif TRYBIAS==2
1271            increaseInThis += oldValue*range;
1272#else
1273            increaseInThis += (oldValue-dualTolerance_)*range;
1274#endif
1275            // goes on swapped list (also means candidates if too many)
1276            spare2[--numberPossiblySwapped]=alpha;
1277            index2[numberPossiblySwapped]=iSequence;
1278            if (fabs(alpha)>bestPivot) {
1279              bestPivot=fabs(alpha);
1280              bestSequence=numberPossiblySwapped;
1281            }
1282          } else {
1283            value = oldValue-upperTheta*alpha;
1284            if (value<-newTolerance && alpha>=acceptablePivot) 
1285              upperTheta = (oldValue+newTolerance)/alpha;
1286            spare2[numberRemaining]=alpha;
1287            index2[numberRemaining++]=iSequence;
1288          }
1289        }
1290      }
1291      swapped[1-iFlip]=numberPossiblySwapped;
1292      interesting[1-iFlip]=numberRemaining;
1293      marker[1-iFlip][0]= max(marker[1-iFlip][0],numberRemaining);
1294      marker[1-iFlip][1]= min(marker[1-iFlip][1],numberPossiblySwapped);
1295     
1296      if (totalThru+thruThis>=fabs(dualOut_)||
1297          increaseInObjective+increaseInThis<0.0) {
1298        // We should be pivoting in this batch
1299        // so compress down to this lot
1300        numberRemaining=0;
1301        for (i=numberColumns_-1;i>=swapped[1-iFlip];i--) {
1302          spare[numberRemaining]=spare2[i];
1303          index[numberRemaining++]=index2[i];
1304        }
1305        interesting[iFlip]=numberRemaining;
1306        int iTry;
1307#define MAXTRY 100
1308        // first get ratio with tolerance
1309        for (iTry=0;iTry<MAXTRY;iTry++) {
1310         
1311          upperTheta=1.0e50;
1312          numberPossiblySwapped = numberColumns_;
1313          numberRemaining = 0;
1314
1315          increaseInThis=0.0; //objective increase in this loop
1316
1317          thruThis=0.0;
1318
1319          spare = array[iFlip];
1320          index = indices[iFlip];
1321          spare2 = array[1-iFlip];
1322          index2 = indices[1-iFlip];
1323     
1324          for (i=0;i<interesting[iFlip];i++) {
1325            int iSequence=index[i];
1326            double alpha=spare[i];
1327            double oldValue = dj_[iSequence];
1328            double value = oldValue-upperTheta*alpha;
1329           
1330            if (alpha < 0.0) {
1331              //at upper bound
1332              if (value>newTolerance) {
1333                if (-alpha>=acceptablePivot) {
1334                  upperTheta = (oldValue-newTolerance)/alpha;
1335                }
1336              }
1337            } else {
1338              // at lower bound
1339              if (value<-newTolerance) {
1340                if (alpha>=acceptablePivot) {
1341                  upperTheta = (oldValue+newTolerance)/alpha;
1342                }
1343              }
1344            }
1345          }
1346          bestPivot=acceptablePivot;
1347          sequenceIn_=-1;
1348          // now choose largest and sum all ones which will go through
1349#define MINIMUMTHETA 1.0e-12
1350          for (i=0;i<interesting[iFlip];i++) {
1351            int iSequence=index[i];
1352            double alpha=spare[i];
1353            double value = dj_[iSequence]-upperTheta*alpha;
1354            double badDj=0.0;
1355           
1356            bool addToSwapped=false;
1357           
1358            if (alpha < 0.0) {
1359              //at upper bound
1360              if (value>=0.0) { 
1361                addToSwapped=true;
1362#if TRYBIAS==1
1363                badDj = -max(dj_[iSequence],0.0);
1364#elif TRYBIAS==2
1365                badDj = -dj_[iSequence];
1366#else
1367                badDj = -dj_[iSequence]-dualTolerance_;
1368#endif
1369              }
1370            } else {
1371              // at lower bound
1372              if (value<=0.0) {
1373                addToSwapped=true;
1374#if TRYBIAS==1
1375                badDj = min(dj_[iSequence],0.0);
1376#elif TRYBIAS==2
1377                badDj = dj_[iSequence];
1378#else
1379                badDj = dj_[iSequence]-dualTolerance_;
1380#endif
1381              }
1382            }
1383            if (!addToSwapped) {
1384              // add to list of remaining
1385              spare2[numberRemaining]=alpha;
1386              index2[numberRemaining++]=iSequence;
1387            } else {
1388              // add to list of swapped
1389              spare2[--numberPossiblySwapped]=alpha;
1390              index2[numberPossiblySwapped]=iSequence;
1391              // select if largest pivot
1392              if (fabs(alpha)>bestPivot) {
1393                sequenceIn_ = numberPossiblySwapped;
1394                bestPivot =  fabs(alpha);
1395                theta_ = dj_[iSequence]/alpha;
1396              }
1397              double range = upper[iSequence] - lower[iSequence];
1398              thruThis += range*fabs(alpha);
1399              increaseInThis += badDj*range;
1400            }
1401          }
1402          swapped[1-iFlip]=numberPossiblySwapped;
1403          interesting[1-iFlip]=numberRemaining;
1404          marker[1-iFlip][0]= max(marker[1-iFlip][0],numberRemaining);
1405          marker[1-iFlip][1]= min(marker[1-iFlip][1],numberPossiblySwapped);
1406          // If we stop now this will be increase in objective (I think)
1407          double increase = (fabs(dualOut_)-totalThru)*theta_;
1408          increase += increaseInObjective;
1409          if (theta_<0.0)
1410            thruThis += fabs(dualOut_); // force using this one
1411          if (increaseInObjective<0.0&&increase<0.0&&lastSequence>=0) {
1412            // back
1413            // We may need to be more careful - we could do by
1414            // switch so we always do fine grained?
1415            bestPivot=0.0;
1416          } else {
1417            // add in
1418            totalThru += thruThis;
1419            increaseInObjective += increaseInThis;
1420          }
1421          if (bestPivot<0.1*bestEverPivot&&
1422              bestEverPivot>1.0e-6&&bestPivot<1.0e-3) {
1423            // back to previous one
1424            sequenceIn_=lastSequence;
1425            // swap regions
1426            iFlip = 1-iFlip;
1427            break;
1428          } else if (sequenceIn_==-1&&upperTheta>largeValue_) {
1429            if (lastPivot>acceptablePivot) {
1430              // back to previous one
1431              sequenceIn_=lastSequence;
1432              // swap regions
1433              iFlip = 1-iFlip;
1434            } else {
1435              // can only get here if all pivots too small
1436            }
1437            break;
1438          } else if (totalThru>=fabs(dualOut_)) {
1439            modifyCosts=true; // fine grain - we can modify costs
1440            break; // no point trying another loop
1441          } else {
1442            lastSequence=sequenceIn_;
1443            if (bestPivot>bestEverPivot)
1444              bestEverPivot=bestPivot;
1445            iFlip = 1 -iFlip;
1446            modifyCosts=true; // fine grain - we can modify costs
1447        }
1448        }
1449        if (iTry==MAXTRY)
1450          iFlip = 1-iFlip; // flip back
1451        break;
1452      } else {
1453        // skip this lot
1454        if (bestPivot>1.0e-3||bestPivot>bestEverPivot) {
1455          bestEverPivot=bestPivot;
1456          lastSequence=bestSequence;
1457        } else {
1458          // keep old swapped
1459          memcpy(array[1-iFlip]+swapped[iFlip],
1460                 array[iFlip]+swapped[iFlip],
1461                 (numberColumns_-swapped[iFlip])*sizeof(double));
1462          memcpy(indices[1-iFlip]+swapped[iFlip],
1463                 indices[iFlip]+swapped[iFlip],
1464                 (numberColumns_-swapped[iFlip])*sizeof(int));
1465          marker[1-iFlip][1] = min(marker[1-iFlip][1],swapped[iFlip]);
1466          swapped[1-iFlip]=swapped[iFlip];
1467        }
1468        increaseInObjective += increaseInThis;
1469        iFlip = 1 - iFlip; // swap regions
1470        tentativeTheta = 2.0*upperTheta;
1471        totalThru += thruThis;
1472      }
1473    }
1474   
1475    // can get here without sequenceIn_ set but with lastSequence
1476    if (sequenceIn_<0&&lastSequence>=0) {
1477      // back to previous one
1478      sequenceIn_=lastSequence;
1479      // swap regions
1480      iFlip = 1-iFlip;
1481    }
1482   
1483    if (sequenceIn_>=0) {
1484      // at this stage sequenceIn_ is just pointer into index array
1485      // flip just so we can use iFlip
1486      iFlip = 1 -iFlip;
1487      spare = array[iFlip];
1488      index = indices[iFlip];
1489      double oldValue;
1490      double alpha = spare[sequenceIn_];
1491      sequenceIn_ = indices[iFlip][sequenceIn_];
1492      oldValue = dj_[sequenceIn_];
1493      theta_ = oldValue/alpha;
1494      if (theta_<MINIMUMTHETA) {
1495        // can't pivot to zero
1496        if (oldValue-MINIMUMTHETA*alpha>=-dualTolerance_) {
1497          theta_=MINIMUMTHETA;
1498        } else if (oldValue-MINIMUMTHETA*alpha>=-newTolerance) {
1499          theta_=MINIMUMTHETA;
1500          thisIncrease=true;
1501        } else {
1502          theta_=(oldValue+newTolerance)/alpha;
1503          assert(theta_>=0.0);
1504          thisIncrease=true;
1505        } 
1506      }
1507      // may need to adjust costs so all dual feasible AND pivoted is exactly 0
1508      if (modifyCosts) {
1509        int i;
1510        double * workRow = rowArray->denseVector();
1511        double * workColumn = columnArray->denseVector();
1512        for (i=numberColumns_-1;i>=swapped[iFlip];i--) {
1513          int iSequence=index[i];
1514          double alpha;
1515          if (iSequence>=numberColumns_)
1516            alpha=workRow[iSequence-numberColumns_];
1517          else
1518            alpha=workColumn[iSequence];
1519          double value = dj_[iSequence]-theta_*alpha;
1520           
1521          // can't be free here
1522         
1523          if (alpha < 0.0) {
1524            //at upper bound
1525            if (value>dualTolerance_) {
1526              thisIncrease=true;
1527#define MODIFYCOST 2
1528#if MODIFYCOST
1529              // modify cost to hit new tolerance
1530              double modification = alpha*theta_-dj_[iSequence]
1531                +newTolerance;
1532              //modification = min(modification,dualTolerance_);
1533              //assert (fabs(modification)<1.0e-7);
1534              dj_[iSequence] += modification;
1535              cost_[iSequence] +=  modification;
1536#endif
1537            }
1538          } else {
1539            // at lower bound
1540            if (-value>dualTolerance_) {
1541              thisIncrease=true;
1542#if MODIFYCOST
1543              // modify cost to hit new tolerance
1544              double modification = alpha*theta_-dj_[iSequence]
1545                -newTolerance;
1546              //modification = max(modification,-dualTolerance_);
1547              //assert (fabs(modification)<1.0e-7);
1548              dj_[iSequence] += modification;
1549              cost_[iSequence] +=  modification;
1550#endif
1551            }
1552          }
1553        }
1554      }
1555    }
1556  }
1557
1558  if (sequenceIn_>=0) {
1559    if (sequenceIn_>=numberColumns_) {
1560      //slack
1561      alpha_ = rowArray->denseVector()[sequenceIn_-numberColumns_];
1562    } else {
1563      // column
1564      alpha_ = columnArray->denseVector()[sequenceIn_];
1565    }
1566    lowerIn_ = lower_[sequenceIn_];
1567    upperIn_ = upper_[sequenceIn_];
1568    valueIn_ = solution_[sequenceIn_];
1569    dualIn_ = dj_[sequenceIn_];
1570
1571    if (numberTimesOptimal_) {
1572      // can we adjust cost back closer to original
1573      //*** add coding
1574    }
1575#if MODIFYCOST>1   
1576    // modify cost to hit zero exactly
1577    // so (dualIn_+modification)==theta_*alpha_
1578    double modification = theta_*alpha_-dualIn_;
1579    dualIn_ += modification;
1580    dj_[sequenceIn_]=dualIn_;
1581    cost_[sequenceIn_] += modification;
1582    //assert (fabs(modification)<1.0e-6);
1583#ifdef CLP_DEBUG
1584    if ((handler_->logLevel()&32)&&fabs(modification)>1.0e-15)
1585      printf("exact %d new cost %g, change %g\n",sequenceIn_,
1586             cost_[sequenceIn_],modification);
1587#endif
1588#endif
1589   
1590    if (alpha_<0.0) {
1591      // as if from upper bound
1592      directionIn_=-1;
1593      upperIn_=valueIn_;
1594    } else {
1595      // as if from lower bound
1596      directionIn_=1;
1597      lowerIn_=valueIn_;
1598    }
1599  }
1600  if (thisIncrease) {
1601    newTolerance = dualTolerance_+1.0e-4*dblParam_[OsiDualTolerance];
1602  }
1603
1604  // clear arrays
1605
1606  for (i=0;i<2;i++) {
1607    memset(array[i],0,marker[i][0]*sizeof(double));
1608    memset(array[i]+marker[i][1],0,
1609           (numberColumns_-marker[i][1])*sizeof(double));
1610  }
1611}
1612/* Checks if tentative optimal actually means unbounded
1613   Returns -3 if not, 2 if is unbounded */
1614int 
1615ClpSimplexDual::checkUnbounded(OsiIndexedVector * ray,
1616                                   OsiIndexedVector * spare,
1617                                   double changeCost)
1618{
1619  int status=2; // say unbounded
1620  factorization_->updateColumn(spare,ray);
1621  // get reduced cost
1622  int i;
1623  int number=ray->getNumElements();
1624  int * index = ray->getIndices();
1625  double * array = ray->denseVector();
1626  for (i=0;i<number;i++) {
1627    int iRow=index[i];
1628    int iPivot=pivotVariable_[iRow];
1629    changeCost -= cost(iPivot)*array[iRow];
1630  }
1631  double way;
1632  if (changeCost>0.0) {
1633    //try going down
1634    way=1.0;
1635  } else if (changeCost<0.0) {
1636    //try going up
1637    way=-1.0;
1638  } else {
1639#ifdef CLP_DEBUG
1640    printf("can't decide on up or down\n");
1641#endif
1642    way=0.0;
1643    status=-3;
1644  }
1645  double movement=1.0e10*way; // some largish number
1646  double zeroTolerance = 1.0e-14*dualBound_;
1647  for (i=0;i<number;i++) {
1648    int iRow=index[i];
1649    int iPivot=pivotVariable_[iRow];
1650    double arrayValue = array[iRow];
1651    if (fabs(arrayValue)<zeroTolerance)
1652      arrayValue=0.0;
1653    double newValue=solution(iPivot)+movement*arrayValue;
1654    if (newValue>upper(iPivot)+primalTolerance_||
1655        newValue<lower(iPivot)-primalTolerance_)
1656      status=-3; // not unbounded
1657  }
1658  if (status==2) {
1659    // create ray
1660    delete [] ray_;
1661    ray_ = new double [numberColumns_];
1662    CoinFillN(ray_,numberColumns_,0.0);
1663    for (i=0;i<number;i++) {
1664      int iRow=index[i];
1665      int iPivot=pivotVariable_[iRow];
1666      double arrayValue = array[iRow];
1667      if (iPivot<numberColumns_&&fabs(arrayValue)>=zeroTolerance)
1668        ray_[iPivot] = way* array[iRow];
1669    }
1670  }
1671  ray->clear();
1672  return status;
1673}
1674/* Checks if finished.  Updates status */
1675void 
1676ClpSimplexDual::statusOfProblemInDual(int & lastCleaned,int type)
1677{
1678  if (type==2) {
1679    // trouble - restore solution
1680    memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
1681    memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
1682           numberRows_*sizeof(double));
1683    memcpy(columnActivityWork_,savedSolution_ ,
1684           numberColumns_*sizeof(double));
1685    forceFactorization_=1; // a bit drastic but ..
1686    changeMade_++; // say something changed
1687  }
1688  int tentativeStatus = problemStatus_;
1689  double changeCost;
1690
1691  if (problemStatus_>-3) {
1692    // factorize
1693    // later on we will need to recover from singularities
1694    // also we could skip if first time
1695    // save dual weights
1696    dualRowPivot_->saveWeights(this,1);
1697    // is factorization okay?
1698    if (internalFactorize(1)) {
1699      // no - restore previous basis
1700      assert (type==1);
1701      changeMade_++; // say something changed
1702      memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
1703      memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
1704             numberRows_*sizeof(double));
1705      memcpy(columnActivityWork_,savedSolution_ ,
1706             numberColumns_*sizeof(double));
1707      // get correct bounds on all variables
1708      double dummyChangeCost=0.0;
1709      changeBounds(true,rowArray_[2],dummyChangeCost);
1710      // throw away change
1711      rowArray_[2]->clear();
1712      forceFactorization_=1; // a bit drastic but ..
1713      type = 2;
1714      assert (internalFactorize(1)==0);
1715    }
1716    problemStatus_=-3;
1717  }
1718  // at this stage status is -3 or -4 if looks infeasible
1719  // get primal and dual solutions
1720  gutsOfSolution(rowActivityWork_,columnActivityWork_);
1721#ifdef CLP_DEBUG
1722  if (!rowScale_&&(handler_->logLevel()&32)) {
1723    double * objectiveSimplex
1724      = copyOfArray(objective_,numberColumns_,0.0);
1725    double * rowObjectiveSimplex
1726      = copyOfArray(rowObjective_,numberRows_,0.0);
1727    int i;
1728    double largest;
1729    largest=0.0;
1730    for (i=0;i<numberRows_;i++) {
1731      rowObjectiveSimplex[i] *= optimizationDirection_;
1732      double difference = fabs(rowObjectiveWork_[i]-rowObjectiveSimplex[i]);
1733      if (difference>largest)
1734        largest=difference;
1735    }
1736    for (i=0;i<numberColumns_;i++) {
1737      objectiveSimplex[i] *= optimizationDirection_;
1738      double difference = fabs(objectiveWork_[i]-objectiveSimplex[i]);
1739      if (difference>largest)
1740        largest=difference;
1741    }
1742    if ((handler_->logLevel()&16))
1743      printf("difference in obj %g\n",largest);
1744    delete [] objectiveSimplex;
1745    delete [] rowObjectiveSimplex;
1746  }
1747#endif
1748  handler_->message(CLP_SIMPLEX_STATUS,messages_)
1749    <<numberIterations_<<objectiveValue();
1750  handler_->printing(sumPrimalInfeasibilities_>0.0)
1751    <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
1752  handler_->printing(sumDualInfeasibilities_>0.0)
1753    <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
1754  handler_->printing(numberDualInfeasibilitiesWithoutFree_
1755                     <numberDualInfeasibilities_)
1756                       <<numberDualInfeasibilities_-
1757    numberDualInfeasibilitiesWithoutFree_;
1758  handler_->message()<<OsiMessageEol;
1759  while (problemStatus_<=-3) {
1760    bool cleanDuals=false;
1761    int numberChangedBounds=0;
1762    int doOriginalTolerance=0;
1763    if ( lastCleaned==numberIterations_)
1764      doOriginalTolerance=1;
1765    // check optimal
1766    if (dualFeasible()||problemStatus_==-4) {
1767      if (primalFeasible()) {
1768        // may be optimal - or may be bounds are wrong
1769        handler_->message(CLP_DUAL_BOUNDS,messages_)
1770          <<dualBound_
1771          <<OsiMessageEol;
1772        // save solution in case unbounded
1773        CoinDisjointCopyN(columnActivityWork_,numberColumns_,
1774                          columnArray_[0]->denseVector());
1775        CoinDisjointCopyN(rowActivityWork_,numberRows_,
1776                          rowArray_[2]->denseVector());
1777        numberChangedBounds=changeBounds(false,rowArray_[0],changeCost);
1778        if (numberChangedBounds<=0) {
1779          //looks optimal - do we need to reset tolerance
1780          if (lastCleaned<numberIterations_&&numberTimesOptimal_<4) {
1781            doOriginalTolerance=2;
1782            numberTimesOptimal_++;
1783            changeMade_++; // say something changed
1784            if (numberTimesOptimal_==1) {
1785              dualTolerance_ = min(dualTolerance_,1.0e-8);
1786              // better to have small tolerance even if slower
1787              factorization_->zeroTolerance(1.0e-15);
1788            }
1789          } else {
1790            problemStatus_=0; // optimal
1791            if (lastCleaned<numberIterations_) {
1792              handler_->message(CLP_SIMPLEX_GIVINGUP,messages_)
1793                <<OsiMessageEol;
1794            }
1795          }
1796        } else {
1797          cleanDuals=true;
1798          if (doOriginalTolerance==1) {
1799            // check unbounded
1800            problemStatus_ = checkUnbounded(rowArray_[0],rowArray_[1],
1801                                            changeCost);
1802            if (problemStatus_==2) {
1803              // it is unbounded - restore solution
1804              // but first add in changes to non-basic
1805              int iColumn;
1806              double * original = columnArray_[0]->denseVector();
1807              for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1808                if(getColumnStatus(iColumn)!= ClpSimplex::basic)
1809                  ray_[iColumn] += 
1810                    columnActivityWork_[iColumn]-original[iColumn];
1811                columnActivityWork_[iColumn] = original[iColumn];
1812              }
1813              CoinDisjointCopyN(rowArray_[2]->denseVector(),numberRows_,
1814                                rowActivityWork_);
1815            }
1816          } else {
1817            doOriginalTolerance=2;
1818            rowArray_[0]->clear();
1819          }
1820        }
1821        CoinFillN(columnArray_[0]->denseVector(),numberColumns_,0.0);
1822        CoinFillN(rowArray_[2]->denseVector(),numberRows_,0.0);
1823      }
1824      if (problemStatus_==-4) {
1825        // may be infeasible - or may be bounds are wrong
1826        handler_->message(CLP_DUAL_CHECKB,messages_)
1827          <<dualBound_
1828          <<OsiMessageEol;
1829        numberChangedBounds=changeBounds(false,NULL,changeCost);
1830        if (numberChangedBounds<=0||dualBound_>1.0e20||
1831            (largestPrimalError_>1.0&&dualBound_>1.0e17)) {
1832          problemStatus_=1; // infeasible
1833        } else {
1834          problemStatus_=-1; //iterate
1835          cleanDuals=true;
1836          doOriginalTolerance=2;
1837          // and delete ray which has been created
1838          delete [] ray_;
1839          ray_ = NULL;
1840        }
1841      }
1842    } else {
1843      cleanDuals=true;
1844    }
1845    if (problemStatus_<0) {
1846      if (doOriginalTolerance==2) {
1847        // put back original tolerance
1848        lastCleaned=numberIterations_;
1849        handler_->message(CLP_DUAL_ORIGINAL,messages_)
1850          <<OsiMessageEol;
1851
1852        perturbation_=102; // stop any perturbations
1853        createRim(4);
1854        // make sure duals are current
1855        computeDuals();
1856        // put back bounds as they were if was optimal
1857        if (doOriginalTolerance==2) {
1858          changeMade_++; // say something changed
1859          changeBounds(true,NULL,changeCost);
1860          cleanDuals=true;
1861        }
1862      }
1863      if (cleanDuals) {
1864        // make sure dual feasible
1865        // look at all rows and columns
1866        rowArray_[0]->clear();
1867        CoinIotaN(rowArray_[0]->getIndices(),numberRows_,0);
1868        rowArray_[0]->setNumElements(numberRows_);
1869        columnArray_[0]->clear();
1870        CoinIotaN(columnArray_[0]->getIndices(),numberColumns_,0);
1871        columnArray_[0]->setNumElements(numberColumns_);
1872        double objectiveChange=0.0;
1873        updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[1],
1874                          0.0,objectiveChange);
1875        // for now - recompute all
1876        gutsOfSolution(rowActivityWork_,columnActivityWork_);
1877        assert(numberDualInfeasibilitiesWithoutFree_==0);
1878        if (numberDualInfeasibilities_) {
1879          // bad free variables
1880          if (primalFeasible()) {
1881            std::cerr<<"Free variable problem?"<<std::endl;
1882            abort(); // what now
1883          }
1884          problemStatus_=-1; // carry on as normal
1885        }
1886      } else {
1887        // iterate
1888        problemStatus_=-1;
1889      }
1890    }
1891  }
1892  if (type==0||type==1) {
1893    if (!type) {
1894      // create save arrays
1895      delete [] saveStatus_;
1896      delete [] savedSolution_;
1897      saveStatus_ = new unsigned char [numberRows_+numberColumns_];
1898      savedSolution_ = new double [numberRows_+numberColumns_];
1899    }
1900    // save arrays
1901    memcpy(saveStatus_,status_,(numberColumns_+numberRows_)*sizeof(char));
1902    memcpy(savedSolution_+numberColumns_ ,rowActivityWork_,
1903           numberRows_*sizeof(double));
1904    memcpy(savedSolution_ ,columnActivityWork_,numberColumns_*sizeof(double));
1905  }
1906
1907  // restore weights (if saved) - also recompute infeasibility list
1908  if (tentativeStatus>-3) 
1909    dualRowPivot_->saveWeights(this,(type <2) ? 2 : 4);
1910  else
1911    dualRowPivot_->saveWeights(this,3);
1912  // unflag all variables (we may want to wait a bit?)
1913  int iRow;
1914  for (iRow=0;iRow<numberRows_;iRow++) {
1915    int iPivot=pivotVariable_[iRow];
1916    clearFlagged(iPivot);
1917  }
1918  if (problemStatus_<0&&!changeMade_) {
1919    problemStatus_=4; // unknown
1920  }
1921
1922}
1923/* While updateDualsInDual sees what effect is of flip
1924   this does actuall flipping.
1925   If change >0.0 then value in array >0.0 => from lower to upper
1926*/
1927void 
1928ClpSimplexDual::flipBounds(OsiIndexedVector * rowArray,
1929                  OsiIndexedVector * columnArray,
1930                  double change)
1931{
1932  double * work;
1933  int number;
1934  int * which;
1935 
1936  int iSection;
1937
1938  for (iSection=0;iSection<2;iSection++) {
1939    int i;
1940    double * solution = solutionRegion(iSection);
1941    double * lower = lowerRegion(iSection);
1942    double * upper = upperRegion(iSection);
1943    int addSequence;
1944    if (!iSection) {
1945      work = rowArray->denseVector();
1946      number = rowArray->getNumElements();
1947      which = rowArray->getIndices();
1948      addSequence = numberColumns_;
1949    } else {
1950      work = columnArray->denseVector();
1951      number = columnArray->getNumElements();
1952      which = columnArray->getIndices();
1953      addSequence = 0;
1954    }
1955   
1956    for (i=0;i<number;i++) {
1957      int iSequence = which[i];
1958#ifndef NDEBUG
1959      double value = work[iSequence]*change;
1960#endif
1961      work[iSequence]=0.0;
1962      Status status = getStatus(iSequence+addSequence);
1963       
1964      switch(status) {
1965       
1966      case ClpSimplex::basic:
1967      case ClpSimplex::isFree:
1968      case ClpSimplex::superBasic:
1969        break;
1970      case ClpSimplex::atUpperBound:
1971        assert (value<=0.0);
1972        // to lower bound
1973        setStatus(iSequence+addSequence,ClpSimplex::atLowerBound);
1974        solution[iSequence] = lower[iSequence];
1975        break;
1976      case ClpSimplex::atLowerBound:
1977        assert (value>=0.0);
1978        // to upper bound
1979        setStatus(iSequence+addSequence,ClpSimplex::atUpperBound);
1980        solution[iSequence] = upper[iSequence];
1981        break;
1982      }
1983    }
1984  }
1985  rowArray->setNumElements(0);
1986  columnArray->setNumElements(0);
1987}
1988// Restores bound to original bound
1989void 
1990ClpSimplexDual::originalBound( int iSequence)
1991{
1992  if (iSequence>=numberColumns_) {
1993    // rows
1994    int iRow = iSequence-numberColumns_;
1995    rowLowerWork_[iRow]=rowLower_[iRow];
1996    rowUpperWork_[iRow]=rowUpper_[iRow];
1997    if (rowScale_) {
1998      if (rowLowerWork_[iRow]>-1.0e50)
1999        rowLowerWork_[iRow] *= rowScale_[iRow];
2000      if (rowUpperWork_[iRow]<1.0e50)
2001        rowUpperWork_[iRow] *= rowScale_[iRow];
2002    }
2003  } else {
2004    // columns
2005    columnLowerWork_[iSequence]=columnLower_[iSequence];
2006    columnUpperWork_[iSequence]=columnUpper_[iSequence];
2007    if (rowScale_) {
2008      double multiplier = 1.0/columnScale_[iSequence];
2009      if (columnLowerWork_[iSequence]>-1.0e50)
2010        columnLowerWork_[iSequence] *= multiplier;
2011      if (columnUpperWork_[iSequence]<1.0e50)
2012        columnUpperWork_[iSequence] *= multiplier;
2013    }
2014  }
2015  setFakeBound(iSequence,ClpSimplex::noFake);
2016}
2017/* As changeBounds but just changes new bounds for a single variable.
2018   Returns true if change */
2019bool 
2020ClpSimplexDual::changeBound( int iSequence)
2021{
2022  // old values
2023  double oldLower=lower_[iSequence];
2024  double oldUpper=upper_[iSequence];
2025  double value=solution_[iSequence];
2026  bool modified=false;
2027  originalBound(iSequence);
2028  // original values
2029  double lowerValue=lower_[iSequence];
2030  double upperValue=upper_[iSequence];
2031  // back to altered values
2032  lower_[iSequence] = oldLower;
2033  upper_[iSequence] = oldUpper;
2034  if (value==oldLower) {
2035    if (upperValue > oldLower + dualBound_) {
2036      upper_[iSequence]=oldLower+dualBound_;
2037      setFakeBound(iSequence,ClpSimplex::upperFake);
2038      modified=true;
2039    }
2040  } else if (value==oldUpper) {
2041    if (lowerValue < oldUpper - dualBound_) {
2042      lower_[iSequence]=oldUpper-dualBound_;
2043      setFakeBound(iSequence,ClpSimplex::lowerFake);
2044      modified=true;
2045    }
2046  } else {
2047    assert(value==oldLower||value==oldUpper);
2048  }
2049  return modified;
2050}
2051// Perturbs problem
2052void 
2053ClpSimplexDual::perturb()
2054{
2055  if (perturbation_>100)
2056    return; //perturbed already
2057  int iRow,iColumn;
2058  // dual perturbation
2059  double perturbation=1.0e-20;
2060  // maximum fraction of cost to perturb
2061  double maximumFraction = 1.0e-4;
2062  if (perturbation_==100) {
2063    perturbation = 1.0e-4;
2064    for (iRow=0;iRow<numberRows_;iRow++) {
2065      double value = fabs(rowActivityWork_[iRow]*rowObjectiveWork_[iRow]);
2066      perturbation = max(perturbation,value);
2067    }
2068    for (iColumn=0;iColumn<numberColumns_;iColumn++) { 
2069      double value = 
2070        fabs(columnActivityWork_[iColumn]*objectiveWork_[iColumn]);
2071      perturbation = max(perturbation,value);
2072    }
2073    perturbation *= 1.0e-8;
2074  } else if (perturbation_<100) {
2075    perturbation = pow(10.0,perturbation_);
2076    // user is in charge
2077    maximumFraction = 1.0e100;
2078  }
2079  // modify costs
2080  handler_->message(CLP_SIMPLEX_PERTURB,messages_)
2081    <<perturbation
2082    <<OsiMessageEol;
2083  for (iRow=0;iRow<numberRows_;iRow++) {
2084    double value = perturbation;
2085    double currentValue = rowObjectiveWork_[iRow];
2086    value = min(value,maximumFraction*fabs(currentValue)+1.0e-6);
2087    if (rowLowerWork_[iRow]>-largeValue_) {
2088      if (fabs(rowLowerWork_[iRow])<fabs(rowUpperWork_[iRow])) 
2089        value *= drand48();
2090      else
2091        value *= -drand48();
2092    } else if (rowUpperWork_[iRow]<largeValue_) {
2093      value *= -drand48();
2094    } else {
2095      value=0.0;
2096    }
2097    rowObjectiveWork_[iRow] += value;
2098  }
2099  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
2100    double value = perturbation;
2101    double currentValue = objectiveWork_[iColumn];
2102    value = min(value,maximumFraction*fabs(currentValue)+1.0e-6);
2103    if (columnLowerWork_[iColumn]>-largeValue_) {
2104      if (fabs(columnLowerWork_[iColumn])<
2105          fabs(columnUpperWork_[iColumn])) 
2106        value *= drand48();
2107      else
2108        value *= -drand48();
2109    } else if (columnUpperWork_[iColumn]<largeValue_) {
2110      value *= -drand48();
2111    } else {
2112      value=0.0;
2113    }
2114    objectiveWork_[iColumn] += value;
2115  }
2116  // say perturbed
2117  perturbation_=102;
2118
2119}
2120
Note: See TracBrowser for help on using the repository browser.