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

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

Adding Clp to development branch

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