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

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

Breaking out whileIterating

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