# source:branches/devel-1/ClpSimplexDual.cpp@7

Last change on this file since 7 was 7, checked in by forrest, 19 years ago

For infeasible problems

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