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

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

Improving reliability using IBM Burlington problems (Primal)

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