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

Last change on this file since 15 was 15, checked in by forrest, 17 years ago

Hope this works from wincvs

Fix error in values pass

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