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

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

Synchronizing

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 78.4 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 "CoinPackedMatrix.hpp"
100#include "CoinIndexedVector.hpp"
101#include "CoinWarmStartBasis.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_[ClpDualTolerance];
240  primalTolerance_=dblParam_[ClpPrimalTolerance];
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    <<CoinMessageEol;
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
281  numberFake_ =0; // Number of variables at fake bounds
282  changeBounds(true,NULL,objectiveChange);
283
284  problemStatus_ = -1;
285  numberIterations_=0;
286
287  int lastCleaned=0; // last time objective or bounds cleaned up
288
289  // number of times we have declared optimality
290  numberTimesOptimal_=0;
291
292  // Progress indicator
293  ClpSimplexProgress progress(this);
294
295  // This says whether to restore things etc
296  int factorType=0;
297  /*
298    Status of problem:
299    0 - optimal
300    1 - infeasible
301    2 - unbounded
302    -1 - iterating
303    -2 - factorization wanted
304    -3 - redo checking without factorization
305    -4 - looks infeasible
306  */
307  while (problemStatus_<0) {
308    // clear
309    for (iRow=0;iRow<4;iRow++) {
310      rowArray_[iRow]->clear();
311    }   
312   
313    for (iColumn=0;iColumn<2;iColumn++) {
314      columnArray_[iColumn]->clear();
315    }   
316
317    // give matrix (and model costs and bounds a chance to be
318    // refreshed (normally null)
319    matrix_->refresh(this);
320    // If getting nowhere - why not give it a kick
321#if 0
322    // does not seem to work too well - do some more work
323    if (perturbation_<101&&numberIterations_>2*(numberRows_+numberColumns_))
324      perturb();
325#endif
326    // may factorize, checks if problem finished
327    statusOfProblemInDual(lastCleaned,factorType,progress);
328
329    // Say good factorization
330    factorType=1;
331    if (saveSparse) {
332      // use default at present
333      factorization_->sparseThreshold(0);
334      factorization_->goSparse();
335    }
336
337    // Do iterations
338    whileIterating();
339  }
340
341  //assert(!numberFake_||problemStatus_); // all bounds should be okay
342  // at present we are leaving factorization around
343  // maybe we should empty it
344  deleteRim();
345  handler_->message(CLP_SIMPLEX_FINISHED+problemStatus_,messages_)
346    <<objectiveValue()
347    <<CoinMessageEol;
348  // Restore any saved stuff
349  perturbation_ = savePerturbation;
350  factorization_->sparseThreshold(saveSparse);
351  dualBound_ = saveDualBound_;
352  return problemStatus_;
353}
354/* Reasons to come out:
355   -1 iterations etc
356   -2 inaccuracy
357   -3 slight inaccuracy (and done iterations)
358   +0 looks optimal (might be unbounded - but we will investigate)
359   +1 looks infeasible
360   +3 max iterations
361 */
362int
363ClpSimplexDual::whileIterating()
364{
365  // status stays at -1 while iterating, >=0 finished, -2 to invert
366  // status -3 to go to top without an invert
367  int returnCode = -1;
368  while (problemStatus_==-1) {
369#ifdef CLP_DEBUG
370    {
371      int i;
372      for (i=0;i<4;i++) {
373        rowArray_[i]->checkClear();
374      }   
375      for (i=0;i<2;i++) {
376        columnArray_[i]->checkClear();
377      }   
378    }     
379#endif
380#if CLP_DEBUG>2
381    // very expensive
382    if (numberIterations_>0&&numberIterations_<-801) {
383      handler_->setLogLevel(63);
384      double saveValue = objectiveValue_;
385      double * saveRow1 = new double[numberRows_];
386      double * saveRow2 = new double[numberRows_];
387      memcpy(saveRow1,rowReducedCost_,numberRows_*sizeof(double));
388      memcpy(saveRow2,rowActivityWork_,numberRows_*sizeof(double));
389      double * saveColumn1 = new double[numberColumns_];
390      double * saveColumn2 = new double[numberColumns_];
391      memcpy(saveColumn1,reducedCostWork_,numberColumns_*sizeof(double));
392      memcpy(saveColumn2,columnActivityWork_,numberColumns_*sizeof(double));
393      gutsOfSolution(rowActivityWork_,columnActivityWork_);
394      printf("xxx %d old obj %g, recomputed %g, sum dual inf %g\n",
395             numberIterations_,
396             saveValue,objectiveValue_,sumDualInfeasibilities_);
397      if (saveValue>objectiveValue_+1.0e-2)
398        printf("**bad**\n");
399      memcpy(rowReducedCost_,saveRow1,numberRows_*sizeof(double));
400      memcpy(rowActivityWork_,saveRow2,numberRows_*sizeof(double));
401      memcpy(reducedCostWork_,saveColumn1,numberColumns_*sizeof(double));
402      memcpy(columnActivityWork_,saveColumn2,numberColumns_*sizeof(double));
403      delete [] saveRow1;
404      delete [] saveRow2;
405      delete [] saveColumn1;
406      delete [] saveColumn2;
407      objectiveValue_=saveValue;
408    }
409#endif
410#ifdef CLP_DEBUG
411    {
412      int iSequence, number=numberRows_+numberColumns_;
413      for (iSequence=0;iSequence<number;iSequence++) {
414        double lowerValue=lower_[iSequence];
415        double upperValue=upper_[iSequence];
416        double value=solution_[iSequence];
417        if(getStatus(iSequence)!=ClpSimplex::basic) {
418          assert(lowerValue>-1.0e20);
419          assert(upperValue<1.0e20);
420        }
421        switch(getStatus(iSequence)) {
422         
423        case ClpSimplex::basic:
424          break;
425        case ClpSimplex::isFree:
426        case ClpSimplex::superBasic:
427          break;
428        case ClpSimplex::atUpperBound:
429          assert (fabs(value-upperValue)<=primalTolerance_) ;
430          break;
431        case ClpSimplex::atLowerBound:
432          assert (fabs(value-lowerValue)<=primalTolerance_) ;
433          break;
434        }
435      }
436    }
437    if(numberIterations_==debugIteration) {
438      printf("dodgy iteration coming up\n");
439    }
440#endif
441    // choose row to go out
442    dualRow();
443    if (pivotRow_>=0) {
444      // we found a pivot row
445      handler_->message(CLP_SIMPLEX_PIVOTROW,messages_)
446        <<pivotRow_
447        <<CoinMessageEol;
448      // check accuracy of weights
449      dualRowPivot_->checkAccuracy();
450      // get sign for finding row of tableau
451      rowArray_[0]->insert(pivotRow_,directionOut_);
452      factorization_->updateColumnTranspose(rowArray_[1],rowArray_[0]);
453      // put row of tableau in rowArray[0] and columnArray[0]
454      matrix_->transposeTimes(this,-1.0,
455                              rowArray_[0],columnArray_[1],columnArray_[0]);
456      // rowArray has pi equivalent
457      // do ratio test
458      dualColumn(rowArray_[0],columnArray_[0],columnArray_[1],
459                 rowArray_[3]);
460      if (sequenceIn_>=0) {
461        // normal iteration
462        // update the incoming column
463        unpack(rowArray_[1]);
464        factorization_->updateColumn(rowArray_[2],rowArray_[1],true);
465        // and update dual weights (can do in parallel - with extra array)
466        dualRowPivot_->updateWeights(rowArray_[0],rowArray_[2],
467                                     rowArray_[1]);
468        // see if update stable
469        double btranAlpha = -alpha_*directionOut_; // for check
470        alpha_=(*rowArray_[1])[pivotRow_];
471#ifdef CLP_DEBUG
472        if ((handler_->logLevel()&32))
473          printf("btran alpha %g, ftran alpha %g\n",btranAlpha,alpha_);
474#endif
475        if (fabs(btranAlpha)<1.0e-12||fabs(alpha_)<1.0e-12||
476            fabs(btranAlpha-alpha_)>1.0e-7*(1.0+fabs(alpha_))) {
477          handler_->message(CLP_DUAL_CHECK,messages_)
478            <<btranAlpha
479            <<alpha_
480            <<CoinMessageEol;
481          dualRowPivot_->unrollWeights();
482          if (factorization_->pivots()) {
483            problemStatus_=-2; // factorize now
484            rowArray_[0]->clear();
485            rowArray_[1]->clear();
486            columnArray_[0]->clear();
487            returnCode=-2;
488            break;
489          } else {
490            // take on more relaxed criterion
491            if (fabs(btranAlpha)<1.0e-12||fabs(alpha_)<1.0e-12||
492                fabs(btranAlpha-alpha_)>1.0e-4*(1.0+fabs(alpha_))) {
493              // need to reject something
494              char x = isColumn(sequenceOut_) ? 'C' :'R';
495              handler_->message(CLP_SIMPLEX_FLAG,messages_)
496                <<x<<sequenceWithin(sequenceOut_)
497                <<CoinMessageEol;
498              setFlagged(sequenceOut_);
499              lastBadIteration_ = numberIterations_; // say be more cautious
500              rowArray_[0]->clear();
501              rowArray_[1]->clear();
502              columnArray_[0]->clear();
503              continue;
504            }
505          }
506        }
507        // update duals BEFORE replaceColumn so can do updateColumn
508        double objectiveChange=0.0;
509        // do duals first as variables may flip bounds
510        // rowArray_[0] and columnArray_[0] may have flips
511        // so use rowArray_[3] for work array from here on
512        int nswapped = 
513          updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[2],theta_,
514                            objectiveChange);
515        // which will change basic solution
516        if (nswapped) {
517#ifdef CLP_DEBUG
518          if ((handler_->logLevel()&16))
519            printf("old dualOut_ %g, v %g, l %g, u %g - new ",
520                   dualOut_,valueOut_,lowerOut_,upperOut_);
521          double oldOut=dualOut_;
522#endif
523          factorization_->updateColumn(rowArray_[3],rowArray_[2],false);
524          dualRowPivot_->updatePrimalSolution(rowArray_[2],
525                                              1.0,objectiveChange);
526         
527          // recompute dualOut_
528          valueOut_ = solution_[sequenceOut_];
529          if (directionOut_<0) {
530            dualOut_ = valueOut_ - upperOut_;
531          } else {
532            dualOut_ = lowerOut_ - valueOut_;
533          }
534#ifdef CLP_DEBUG
535          if ((handler_->logLevel()&16))
536            printf("%g\n",dualOut_);
537          assert(dualOut_<=oldOut);
538#endif
539          if(dualOut_<0.0&&factorization_->pivots()) {
540            // going backwards - factorize
541            dualRowPivot_->unrollWeights();
542            problemStatus_=-2; // factorize now
543            returnCode=-2;
544            break;
545          }
546        }
547        // amount primal will move
548        double movement = -dualOut_*directionOut_/alpha_;
549        // so objective should increase by fabs(dj)*movement
550        // but we already have objective change - so check will be good
551        if (objectiveChange+fabs(movement*dualIn_)<-1.0e-5) {
552#ifdef CLP_DEBUG
553          if (handler_->logLevel()&32)
554            printf("movement %g, swap change %g, rest %g  * %g\n",
555                   objectiveChange+fabs(movement*dualIn_),
556                   objectiveChange,movement,dualIn_);
557#endif
558          if(factorization_->pivots()>5) {
559            // going backwards - factorize
560            dualRowPivot_->unrollWeights();
561            problemStatus_=-2; // factorize now
562            returnCode=-2;
563            break;
564          }
565        }
566        // if stable replace in basis
567        int updateStatus = factorization_->replaceColumn(rowArray_[2],
568                                                         pivotRow_,
569                                                         alpha_);
570        if (updateStatus==1) {
571          // slight error
572          if (factorization_->pivots()>5) {
573            problemStatus_=-2; // factorize now
574            returnCode=-3;
575          }
576        } else if (updateStatus==2) {
577          // major error
578          dualRowPivot_->unrollWeights();
579          // later we may need to unwind more e.g. fake bounds
580          if (factorization_->pivots()) {
581            problemStatus_=-2; // factorize now
582            returnCode=-2;
583            break;
584          } else {
585            // need to reject something
586            char x = isColumn(sequenceOut_) ? 'C' :'R';
587            handler_->message(CLP_SIMPLEX_FLAG,messages_)
588              <<x<<sequenceWithin(sequenceOut_)
589              <<CoinMessageEol;
590            setFlagged(sequenceOut_);
591            lastBadIteration_ = numberIterations_; // say be more cautious
592            rowArray_[0]->clear();
593            rowArray_[1]->clear();
594            columnArray_[0]->clear();
595            continue;
596          }
597        } else if (updateStatus==3) {
598          // out of memory
599          // increase space if not many iterations
600          if (factorization_->pivots()<
601              0.5*factorization_->maximumPivots()&&
602              factorization_->pivots()<200)
603            factorization_->areaFactor(
604                                       factorization_->areaFactor() * 1.1);
605          problemStatus_=-2; // factorize now
606        } 
607        // update primal solution
608        if (theta_<0.0) {
609#ifdef CLP_DEBUG
610          if (handler_->logLevel()&32)
611            printf("negative theta %g\n",theta_);
612#endif
613          theta_=0.0;
614        }
615        // do actual flips
616        flipBounds(rowArray_[0],columnArray_[0],theta_);
617        dualRowPivot_->updatePrimalSolution(rowArray_[1],
618                                            movement,
619                                            objectiveChange);
620#ifdef CLP_DEBUG
621        double oldobj=objectiveValue_;
622#endif
623        int whatNext=housekeeping(objectiveChange);
624        // and set bounds correctly
625        originalBound(sequenceIn_); 
626        changeBound(sequenceOut_);
627#ifdef CLP_DEBUG
628        if (objectiveValue_<oldobj-1.0e-5&&(handler_->logLevel()&16))
629          printf("obj backwards %g %g\n",objectiveValue_,oldobj);
630#endif
631        if (whatNext==1) {
632          problemStatus_ =-2; // refactorize
633        } else if (whatNext==2) {
634          // maximum iterations or equivalent
635          problemStatus_= 3;
636          returnCode=3;
637          break;
638        }
639      } else {
640        // no incoming column is valid
641#ifdef CLP_DEBUG
642        if (handler_->logLevel()&32)
643          printf("** no column pivot\n");
644#endif
645        if (factorization_->pivots()<5) {
646          problemStatus_=-4; //say looks infeasible
647          // create ray anyway
648          delete [] ray_;
649          ray_ = new double [ numberRows_];
650          CoinDisjointCopyN(rowArray_[0]->denseVector(),numberRows_,ray_);
651        }
652        rowArray_[0]->clear();
653        columnArray_[0]->clear();
654        returnCode=1;
655        break;
656      }
657    } else {
658      // no pivot row
659#ifdef CLP_DEBUG
660      if (handler_->logLevel()&32)
661        printf("** no row pivot\n");
662#endif
663      if (!factorization_->pivots()) {
664        // may have crept through - so may be optimal
665        //problemStatus_=-5; //say looks unbounded
666        problemStatus_=0;
667        // check any flagged variables
668        int iRow;
669        for (iRow=0;iRow<numberRows_;iRow++) {
670          int iPivot=pivotVariable_[iRow];
671          if (flagged(iPivot))
672            break;
673        }
674        if (iRow<numberRows_) {
675#ifdef CLP_DEBUG
676          std::cerr<<"Flagged variables at end - infeasible?"<<std::endl;
677#endif
678          problemStatus_=-4; //say looks infeasible
679          // create ray anyway
680          delete [] ray_;
681          ray_ = new double [ numberRows_];
682          CoinDisjointCopyN(rowArray_[0]->denseVector(),numberRows_,ray_);
683        }
684      }
685      returnCode=0;
686      break;
687    }
688  }
689  return returnCode;
690}
691/* The duals are updated by the given arrays.
692   Returns number of infeasibilities.
693   rowArray and columnarray will have flipped
694   The output vector has movement (row length array) */
695int
696ClpSimplexDual::updateDualsInDual(CoinIndexedVector * rowArray,
697                                  CoinIndexedVector * columnArray,
698                                  CoinIndexedVector * outputArray,
699                                  double theta,
700                                  double & objectiveChange)
701{
702 
703  outputArray->clear();
704 
705  double * work;
706  int number;
707  int * which;
708 
709  int numberInfeasibilities=0;
710  int numberRowInfeasibilities=0;
711 
712  // see whether we will be doing full recompute
713  bool fullRecompute= (rowArray->getNumElements()==numberRows_&&
714                       columnArray->getNumElements()==numberColumns_);
715  int numberAtFake=0;
716 
717  // use a tighter tolerance except for all being okay
718  double tolerance = dualTolerance_;
719 
720  double changeObj=0.0;
721 
722  int iSection;
723 
724  for (iSection=0;iSection<2;iSection++) {
725    int i;
726    double * solution = solutionRegion(iSection);
727    double * reducedCost = djRegion(iSection);
728    double * lower = lowerRegion(iSection);
729    double * upper = upperRegion(iSection);
730    double * cost = costRegion(iSection);
731    int addSequence;
732    if (!iSection) {
733      addSequence = numberColumns_;
734      work = rowArray->denseVector();
735      number = rowArray->getNumElements();
736      which = rowArray->getIndices();
737    } else {
738      // set number of infeasibilities in row array
739      addSequence=0;
740      numberRowInfeasibilities=numberInfeasibilities;
741      rowArray->setNumElements(numberInfeasibilities);
742      numberInfeasibilities=0;
743      work = columnArray->denseVector();
744      number = columnArray->getNumElements();
745      which = columnArray->getIndices();
746    }
747   
748    for (i=0;i<number;i++) {
749      int iSequence = which[i];
750      double alphaI = work[iSequence];
751      double value = reducedCost[iSequence]-theta*alphaI;
752      work[iSequence]=0.0;
753      reducedCost[iSequence]=value;
754     
755      if (!fixed(iSequence+addSequence)) {
756        double movement=0.0;
757        FakeBound bound = getFakeBound(iSequence+addSequence);
758        Status status = getStatus(iSequence+addSequence);
759
760        switch(status) {
761         
762        case ClpSimplex::basic:
763        case ClpSimplex::superBasic:
764          break;
765        case ClpSimplex::isFree:
766          if (fabs(value)>tolerance) { 
767#ifdef CLP_DEBUG
768            if (handler_->logLevel()&32)
769              printf("%d %d, free has dj of %g, alpha %g\n",
770                     iSection,iSequence,value,alphaI);
771#endif
772          }
773          break;
774        case ClpSimplex::atUpperBound:
775          if (value>tolerance) {
776            // to lower bound (if swap)
777            // put back alpha
778            which[numberInfeasibilities++]=iSequence;
779            work[iSequence]=alphaI;
780            movement = lower[iSequence]-upper[iSequence];
781#ifdef CLP_DEBUG
782            if ((handler_->logLevel()&32))
783              printf("%d %d, new dj %g, alpha %g, movement %g\n",
784                   iSection,iSequence,value,alphaI,movement);
785#endif
786            changeObj += movement*cost[iSequence];
787            if (bound==ClpSimplexDual::bothFake||
788                bound==ClpSimplexDual::lowerFake) 
789              numberAtFake++;
790          } else if (fullRecompute) {
791            // at correct bound
792            if (bound==ClpSimplexDual::bothFake||
793                bound==ClpSimplexDual::upperFake) {
794              // but flip if dj would allow
795              if (bound==ClpSimplexDual::upperFake&&
796                  value>=-tolerance) {
797                movement = lower[iSequence]-upper[iSequence];
798                setStatus(iSequence+addSequence,ClpSimplex::atLowerBound);
799                solution[iSequence] = lower[iSequence];
800                changeObj += movement*cost[iSequence];
801              } else {
802                numberAtFake++;
803              }
804            }
805          }
806          break;
807        case ClpSimplex::atLowerBound:
808          if (value<-tolerance) {
809            // to upper bound
810            // put back alpha
811            which[numberInfeasibilities++]=iSequence;
812            work[iSequence]=alphaI;
813            movement = upper[iSequence] - lower[iSequence];
814#ifdef CLP_DEBUG
815            if ((handler_->logLevel()&32))
816              printf("%d %d, new dj %g, alpha %g, movement %g\n",
817                   iSection,iSequence,value,alphaI,movement);
818#endif
819            changeObj += movement*cost[iSequence];
820            if (bound==ClpSimplexDual::bothFake||
821                bound==ClpSimplexDual::upperFake) 
822              numberAtFake++;
823          } else if (fullRecompute) {
824            // at correct bound
825            if (bound==ClpSimplexDual::bothFake||
826                bound==ClpSimplexDual::lowerFake) {
827              // but flip if dj would allow
828              if (bound==ClpSimplexDual::lowerFake&&
829                  value<=tolerance) {
830                movement = upper[iSequence] - lower[iSequence];
831                setStatus(iSequence+addSequence,ClpSimplex::atUpperBound);
832                solution[iSequence] = upper[iSequence];
833                changeObj += movement*cost[iSequence];
834              } else {
835                numberAtFake++;
836              }
837            }
838          }
839          break;
840        }
841        if (!fullRecompute) {
842          if (movement) {
843            if (!iSection) {
844              // row (sign ?)
845              outputArray->quickAdd(iSequence,-movement);
846            } else {
847              matrix_->add(this,outputArray,iSequence,movement);
848            }
849          }
850        }
851      }
852    }
853  }
854#ifdef CLP_DEBUG
855  if (fullRecompute&&numberAtFake&&(handler_->logLevel()&16)!=0) 
856    printf("%d fake after full update\n",numberAtFake);
857#endif
858  outputArray->stopQuickAdd();
859  // set number of infeasibilities
860  columnArray->setNumElements(numberInfeasibilities);
861  numberInfeasibilities += numberRowInfeasibilities;
862  if (fullRecompute) {
863    // do actual flips
864    flipBounds(rowArray,columnArray,theta);
865    numberFake_ = numberAtFake;
866  }
867  objectiveChange += changeObj;
868  return numberInfeasibilities;
869}
870/*
871   Chooses dual pivot row
872   Would be faster with separate region to scan
873   and will have this (with square of infeasibility) when steepest
874   For easy problems we can just choose one of the first rows we look at
875*/
876void
877ClpSimplexDual::dualRow()
878{
879  // get pivot row using whichever method it is
880  pivotRow_=dualRowPivot_->pivotRow();
881  if (pivotRow_>=0) {
882    int iPivot=pivotVariable_[pivotRow_];
883    sequenceOut_ = iPivot;
884    if (iPivot>=numberColumns_) {
885      // slack
886      iPivot-=numberColumns_;
887      valueOut_=rowActivityWork_[iPivot];
888      lowerOut_=rowLowerWork_[iPivot];
889      upperOut_=rowUpperWork_[iPivot];
890    } else {
891      // column
892      valueOut_=columnActivityWork_[iPivot];
893      lowerOut_=columnLowerWork_[iPivot];
894      upperOut_=columnUpperWork_[iPivot];
895    }
896    // if we have problems we could try other way and hope we get a
897    // zero pivot?
898    if (valueOut_>upperOut_) {
899      directionOut_ = -1;
900      dualOut_ = valueOut_ - upperOut_;
901    } else {
902      directionOut_ = 1;
903      dualOut_ = lowerOut_ - valueOut_;
904    }
905#ifdef CLP_DEBUG
906    assert(dualOut_>=0.0);
907#endif
908  }
909  return ;
910}
911// Checks if any fake bounds active - if so returns number and modifies
912// dualBound_ and everything.
913// Free variables will be left as free
914// Returns number of bounds changed if >=0
915// Returns -1 if not initialize and no effect
916// Fills in changeVector which can be used to see if unbounded
917// and cost of change vector
918int
919ClpSimplexDual::changeBounds(bool initialize,
920                                 CoinIndexedVector * outputArray,
921                                 double & changeCost)
922{
923  if (!initialize) {
924    int numberInfeasibilities;
925    double newBound;
926    newBound = 5.0*dualBound_;
927    numberInfeasibilities=0;
928    changeCost=0.0;
929    // put back original bounds and then check
930    createRim(3);
931    int iSequence;
932    // bounds will get bigger - just look at ones at bounds
933    for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
934      double lowerValue=lower_[iSequence];
935      double upperValue=upper_[iSequence];
936      double value=solution_[iSequence];
937      setFakeBound(iSequence,ClpSimplexDual::noFake);
938      switch(getStatus(iSequence)) {
939
940      case ClpSimplex::basic:
941        break;
942      case ClpSimplex::isFree:
943      case ClpSimplex::superBasic:
944        break;
945      case ClpSimplex::atUpperBound:
946        if (fabs(value-upperValue)>primalTolerance_) 
947          numberInfeasibilities++;
948        break;
949      case ClpSimplex::atLowerBound:
950        if (fabs(value-lowerValue)>primalTolerance_) 
951          numberInfeasibilities++;
952        break;
953      }
954    }
955    if (numberInfeasibilities) {
956      int iSequence;
957      for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
958        double lowerValue=lower_[iSequence];
959        double upperValue=upper_[iSequence];
960        double newLowerValue;
961        double newUpperValue;
962        Status status = getStatus(iSequence);
963        if (status==ClpSimplex::atUpperBound||
964            status==ClpSimplex::atLowerBound) {
965          double value = solution_[iSequence];
966          if (value-lowerValue<=upperValue-value) {
967            newLowerValue = max(lowerValue,value-0.666667*newBound);
968            newUpperValue = min(upperValue,newLowerValue+newBound);
969          } else {
970            newUpperValue = min(upperValue,value+0.666667*newBound);
971            newLowerValue = max(lowerValue,newUpperValue-newBound);
972          }
973          lower_[iSequence]=newLowerValue;
974          upper_[iSequence]=newUpperValue;
975          if (newLowerValue > lowerValue) {
976            if (newUpperValue < upperValue) 
977              setFakeBound(iSequence,ClpSimplexDual::bothFake);
978            else
979              setFakeBound(iSequence,ClpSimplexDual::lowerFake);
980          } else {
981            if (newUpperValue < upperValue) 
982              setFakeBound(iSequence,ClpSimplexDual::upperFake);
983          }
984          if (status==ClpSimplex::atUpperBound)
985            solution_[iSequence] = newUpperValue;
986          else
987            solution_[iSequence] = newLowerValue;
988          double movement = solution_[iSequence] - value;
989          if (movement&&outputArray) {
990            if (iSequence>=numberColumns_) {
991              outputArray->quickAdd(iSequence,-movement);
992              changeCost += movement*cost_[iSequence];
993            } else {
994              matrix_->add(this,outputArray,iSequence,movement);
995              changeCost += movement*cost_[iSequence];
996            }
997          }
998        }
999      }
1000      dualBound_ = newBound;
1001      if (outputArray)
1002        outputArray->stopQuickAdd();
1003    } else {
1004      numberInfeasibilities=-1;
1005    }
1006    return numberInfeasibilities;
1007  } else {
1008    int iSequence;
1009     
1010    for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
1011      Status status = getStatus(iSequence);
1012      if (status==ClpSimplex::atUpperBound||
1013          status==ClpSimplex::atLowerBound) {
1014        double lowerValue=lower_[iSequence];
1015        double upperValue=upper_[iSequence];
1016        double value = solution_[iSequence];
1017        if (lowerValue>-largeValue_||upperValue<largeValue_) {
1018          if (lowerValue-value>-0.5*dualBound_||
1019              upperValue-value<0.5*dualBound_) {
1020            if (fabs(lowerValue-value)<=fabs(upperValue-value)) {
1021              if (upperValue > lowerValue + dualBound_) {
1022                upper_[iSequence]=lowerValue+dualBound_;
1023                setFakeBound(iSequence,ClpSimplexDual::upperFake);
1024              }
1025            } else {
1026              if (lowerValue < upperValue - dualBound_) {
1027                lower_[iSequence]=upperValue-dualBound_;
1028                setFakeBound(iSequence,ClpSimplexDual::lowerFake);
1029              }
1030            }
1031          } else {
1032            lower_[iSequence]=-0.5*dualBound_;
1033            upper_[iSequence]= 0.5*dualBound_;
1034            setFakeBound(iSequence,ClpSimplexDual::bothFake);
1035          }
1036        }
1037      }
1038    }
1039    return 1;
1040  }
1041}
1042/*
1043   Row array has row part of pivot row (as duals so sign may be switched)
1044   Column array has column part.
1045   This chooses pivot column.
1046   Spare array will be needed when we start getting clever.
1047   We will check for basic so spare array will never overflow.
1048   If necessary will modify costs
1049*/
1050void
1051ClpSimplexDual::dualColumn(CoinIndexedVector * rowArray,
1052                       CoinIndexedVector * columnArray,
1053                       CoinIndexedVector * spareArray,
1054                       CoinIndexedVector * spareArray2)
1055{
1056  double * work;
1057  int number;
1058  int * which;
1059  double * reducedCost;
1060
1061  int iSection;
1062
1063  sequenceIn_=-1;
1064  int numberPossiblySwapped=0;
1065  int numberRemaining=0;
1066 
1067  double totalThru=0.0; // for when variables flip
1068  double acceptablePivot=1.0e-7;
1069  if (factorization_->pivots())
1070    acceptablePivot=1.0e-5; // if we have iterated be more strict
1071  double bestEverPivot=acceptablePivot;
1072  int lastSequence = -1;
1073  double lastPivot=0.0;
1074  double upperTheta;
1075  double newTolerance = dualTolerance_;
1076  // will we need to increase tolerance
1077  bool thisIncrease=false;
1078  // If we think we need to modify costs (not if something from broad sweep)
1079  bool modifyCosts=false;
1080  // Increase in objective due to swapping bounds (may be negative)
1081  double increaseInObjective=0.0;
1082
1083  // use spareArrays to put ones looked at in
1084  // we are going to flip flop between
1085  int iFlip = 0;
1086  // Possible list of pivots
1087  int interesting[2];
1088  // where possible swapped ones are
1089  int swapped[2];
1090  // for zeroing out arrays after
1091  int marker[2][2];
1092  // pivot elements
1093  double * array[2], * spare, * spare2;
1094  // indices
1095  int * indices[2], * index, * index2;
1096  spareArray->clear();
1097  spareArray2->clear();
1098  array[0] = spareArray->denseVector();
1099  indices[0] = spareArray->getIndices();
1100  spare = array[0];
1101  index = indices[0];
1102  array[1] = spareArray2->denseVector();
1103  indices[1] = spareArray2->getIndices();
1104  int i;
1105  double * lower;
1106  double * upper;
1107
1108  // initialize lists
1109  for (i=0;i<2;i++) {
1110    interesting[i]=0;
1111    swapped[i]=numberColumns_;
1112    marker[i][0]=0;
1113    marker[i][1]=numberColumns_;
1114  }
1115
1116  /*
1117    First we get a list of possible pivots.  We can also see if the
1118    problem looks infeasible or whether we want to pivot in free variable.
1119    This may make objective go backwards but can only happen a finite
1120    number of times and I do want free variables basic.
1121
1122    Then we flip back and forth.  At the start of each iteration
1123    interesting[iFlip] should have possible candidates and swapped[iFlip]
1124    will have pivots if we decide to take a previous pivot.
1125    At end of each iteration interesting[1-iFlip] should have
1126    candidates if we go through this theta and swapped[1-iFlip]
1127    pivots if we don't go through.
1128
1129    At first we increase theta and see what happens.  We start
1130    theta at a reasonable guess.  If in right area then we do bit by bit.
1131
1132   */
1133
1134  // do first pass to get possibles
1135  // We can also see if infeasible or pivoting on free
1136  double tentativeTheta = 1.0e22;
1137  upperTheta = 1.0e31;
1138  double freePivot = acceptablePivot;
1139  for (iSection=0;iSection<2;iSection++) {
1140
1141    int addSequence;
1142
1143    if (!iSection) {
1144      lower = rowLowerWork_;
1145      upper = rowUpperWork_;
1146      work = rowArray->denseVector();
1147      number = rowArray->getNumElements();
1148      which = rowArray->getIndices();
1149      reducedCost = rowReducedCost_;
1150      addSequence = numberColumns_;
1151    } else {
1152      lower = columnLowerWork_;
1153      upper = columnUpperWork_;
1154      work = columnArray->denseVector();
1155      number = columnArray->getNumElements();
1156      which = columnArray->getIndices();
1157      reducedCost = reducedCostWork_;
1158      addSequence = 0;
1159    }
1160   
1161    for (i=0;i<number;i++) {
1162      int iSequence = which[i];
1163      double alpha = work[iSequence];
1164      if (fixed(iSequence+addSequence)||!alpha) 
1165        continue; // skip fixed ones or (zeroed out)
1166      double oldValue = reducedCost[iSequence];
1167      double value = oldValue-tentativeTheta*alpha;
1168      int keep = 0;
1169
1170      switch(getStatus(iSequence+addSequence)) {
1171         
1172      case ClpSimplex::basic:
1173        break;
1174      case ClpSimplex::isFree:
1175      case ClpSimplex::superBasic:
1176        if (oldValue>dualTolerance_) {
1177          if (value<-newTolerance) 
1178            keep = 2;
1179        } else if (oldValue<-dualTolerance_) {
1180          if (value>newTolerance) 
1181            keep = 2;
1182        } else {
1183          if (alpha>=acceptablePivot) 
1184            keep = 2;
1185          else if (-alpha>=acceptablePivot) 
1186            keep = 2;
1187        }
1188        break;
1189      case ClpSimplex::atUpperBound:
1190        assert (oldValue<=dualTolerance_*1.0001);
1191        if (value>newTolerance) {
1192          keep = 1;
1193          value = oldValue-upperTheta*alpha;
1194          if (value>newTolerance && -alpha>=acceptablePivot) 
1195            upperTheta = (oldValue-newTolerance)/alpha;
1196        }
1197        break;
1198      case ClpSimplex::atLowerBound:
1199        assert (oldValue>=-dualTolerance_*1.0001);
1200        if (value<-newTolerance) {
1201          keep = 1;
1202          value = oldValue-upperTheta*alpha;
1203          if (value<-newTolerance && alpha>=acceptablePivot) 
1204            upperTheta = (oldValue+newTolerance)/alpha;
1205        }
1206        break;
1207      }
1208      if (keep) {
1209        if (keep==2) {
1210          // free - choose largest
1211          if (fabs(alpha)>freePivot) {
1212            freePivot=fabs(alpha);
1213            sequenceIn_ = iSequence + addSequence;
1214            theta_=oldValue/alpha;
1215          }
1216        } else {
1217          // add to list
1218          spare[numberRemaining]=alpha;
1219          index[numberRemaining++]=iSequence+addSequence;
1220        }
1221      }
1222    }
1223  }
1224  interesting[0]=numberRemaining;
1225  marker[0][0] = numberRemaining;
1226
1227  if (!numberRemaining)
1228    return; // Looks infeasible
1229
1230  if (sequenceIn_>=0) {
1231    // free variable - always choose
1232  } else {
1233
1234    theta_=1.0e50;
1235    // now flip flop between spare arrays until reasonable theta
1236    tentativeTheta = max(10.0*upperTheta,1.0e-7);
1237   
1238    // loops increasing tentative theta until can't go through
1239   
1240    while (tentativeTheta < 1.0e22) {
1241      double thruThis = 0.0;
1242     
1243      double bestPivot=acceptablePivot;
1244      int bestSequence=-1;
1245     
1246      numberPossiblySwapped = numberColumns_;
1247      numberRemaining = 0;
1248
1249      upperTheta = 1.0e50;
1250
1251      spare = array[iFlip];
1252      index = indices[iFlip];
1253      spare2 = array[1-iFlip];
1254      index2 = indices[1-iFlip];
1255
1256      // try 3 different ways
1257      // 1 bias increase by ones with slightly wrong djs
1258      // 2 bias by all
1259      // 3 bias by all - tolerance (doesn't seem very good)
1260#define TRYBIAS 1
1261
1262
1263      double increaseInThis=0.0; //objective increase in this loop
1264     
1265      for (i=0;i<interesting[iFlip];i++) {
1266        int iSequence = index[i];
1267        double alpha = spare[i];
1268        double oldValue = dj_[iSequence];
1269        double value = oldValue-tentativeTheta*alpha;
1270
1271        if (alpha < 0.0) {
1272          //at upper bound
1273          if (value>newTolerance) {
1274            double range = upper_[iSequence] - lower_[iSequence];
1275            thruThis -= range*alpha;
1276#if TRYBIAS==1
1277            if (oldValue>0.0)
1278              increaseInThis -= oldValue*range;
1279#elif TRYBIAS==2
1280            increaseInThis -= oldValue*range;
1281#else
1282            increaseInThis -= (oldValue+dualTolerance_)*range;
1283#endif
1284            // goes on swapped list (also means candidates if too many)
1285            spare2[--numberPossiblySwapped]=alpha;
1286            index2[numberPossiblySwapped]=iSequence;
1287            if (fabs(alpha)>bestPivot) {
1288              bestPivot=fabs(alpha);
1289              bestSequence=numberPossiblySwapped;
1290            }
1291          } else {
1292            value = oldValue-upperTheta*alpha;
1293            if (value>newTolerance && -alpha>=acceptablePivot) 
1294              upperTheta = (oldValue-newTolerance)/alpha;
1295            spare2[numberRemaining]=alpha;
1296            index2[numberRemaining++]=iSequence;
1297          }
1298        } else {
1299          // at lower bound
1300          if (value<-newTolerance) {
1301            double range = upper_[iSequence] - lower_[iSequence];
1302            thruThis += range*alpha;
1303            //?? is this correct - and should we look at good ones
1304#if TRYBIAS==1
1305            if (oldValue<0.0)
1306              increaseInThis += oldValue*range;
1307#elif TRYBIAS==2
1308            increaseInThis += oldValue*range;
1309#else
1310            increaseInThis += (oldValue-dualTolerance_)*range;
1311#endif
1312            // goes on swapped list (also means candidates if too many)
1313            spare2[--numberPossiblySwapped]=alpha;
1314            index2[numberPossiblySwapped]=iSequence;
1315            if (fabs(alpha)>bestPivot) {
1316              bestPivot=fabs(alpha);
1317              bestSequence=numberPossiblySwapped;
1318            }
1319          } else {
1320            value = oldValue-upperTheta*alpha;
1321            if (value<-newTolerance && alpha>=acceptablePivot) 
1322              upperTheta = (oldValue+newTolerance)/alpha;
1323            spare2[numberRemaining]=alpha;
1324            index2[numberRemaining++]=iSequence;
1325          }
1326        }
1327      }
1328      swapped[1-iFlip]=numberPossiblySwapped;
1329      interesting[1-iFlip]=numberRemaining;
1330      marker[1-iFlip][0]= max(marker[1-iFlip][0],numberRemaining);
1331      marker[1-iFlip][1]= min(marker[1-iFlip][1],numberPossiblySwapped);
1332     
1333      if (totalThru+thruThis>=fabs(dualOut_)||
1334          increaseInObjective+increaseInThis<0.0) {
1335        // We should be pivoting in this batch
1336        // so compress down to this lot
1337        numberRemaining=0;
1338        for (i=numberColumns_-1;i>=swapped[1-iFlip];i--) {
1339          spare[numberRemaining]=spare2[i];
1340          index[numberRemaining++]=index2[i];
1341        }
1342        interesting[iFlip]=numberRemaining;
1343        int iTry;
1344#define MAXTRY 100
1345        // first get ratio with tolerance
1346        for (iTry=0;iTry<MAXTRY;iTry++) {
1347         
1348          upperTheta=1.0e50;
1349          numberPossiblySwapped = numberColumns_;
1350          numberRemaining = 0;
1351
1352          increaseInThis=0.0; //objective increase in this loop
1353
1354          thruThis=0.0;
1355
1356          spare = array[iFlip];
1357          index = indices[iFlip];
1358          spare2 = array[1-iFlip];
1359          index2 = indices[1-iFlip];
1360     
1361          for (i=0;i<interesting[iFlip];i++) {
1362            int iSequence=index[i];
1363            double alpha=spare[i];
1364            double oldValue = dj_[iSequence];
1365            double value = oldValue-upperTheta*alpha;
1366           
1367            if (alpha < 0.0) {
1368              //at upper bound
1369              if (value>newTolerance) {
1370                if (-alpha>=acceptablePivot) {
1371                  upperTheta = (oldValue-newTolerance)/alpha;
1372                }
1373              }
1374            } else {
1375              // at lower bound
1376              if (value<-newTolerance) {
1377                if (alpha>=acceptablePivot) {
1378                  upperTheta = (oldValue+newTolerance)/alpha;
1379                }
1380              }
1381            }
1382          }
1383          bestPivot=acceptablePivot;
1384          sequenceIn_=-1;
1385          // now choose largest and sum all ones which will go through
1386#define MINIMUMTHETA 1.0e-12
1387          for (i=0;i<interesting[iFlip];i++) {
1388            int iSequence=index[i];
1389            double alpha=spare[i];
1390            double value = dj_[iSequence]-upperTheta*alpha;
1391            double badDj=0.0;
1392           
1393            bool addToSwapped=false;
1394           
1395            if (alpha < 0.0) {
1396              //at upper bound
1397              if (value>=0.0) { 
1398                addToSwapped=true;
1399#if TRYBIAS==1
1400                badDj = -max(dj_[iSequence],0.0);
1401#elif TRYBIAS==2
1402                badDj = -dj_[iSequence];
1403#else
1404                badDj = -dj_[iSequence]-dualTolerance_;
1405#endif
1406              }
1407            } else {
1408              // at lower bound
1409              if (value<=0.0) {
1410                addToSwapped=true;
1411#if TRYBIAS==1
1412                badDj = min(dj_[iSequence],0.0);
1413#elif TRYBIAS==2
1414                badDj = dj_[iSequence];
1415#else
1416                badDj = dj_[iSequence]-dualTolerance_;
1417#endif
1418              }
1419            }
1420            if (!addToSwapped) {
1421              // add to list of remaining
1422              spare2[numberRemaining]=alpha;
1423              index2[numberRemaining++]=iSequence;
1424            } else {
1425              // add to list of swapped
1426              spare2[--numberPossiblySwapped]=alpha;
1427              index2[numberPossiblySwapped]=iSequence;
1428              // select if largest pivot
1429              if (fabs(alpha)>bestPivot) {
1430                sequenceIn_ = numberPossiblySwapped;
1431                bestPivot =  fabs(alpha);
1432                theta_ = dj_[iSequence]/alpha;
1433              }
1434              double range = upper[iSequence] - lower[iSequence];
1435              thruThis += range*fabs(alpha);
1436              increaseInThis += badDj*range;
1437            }
1438          }
1439          swapped[1-iFlip]=numberPossiblySwapped;
1440          interesting[1-iFlip]=numberRemaining;
1441          marker[1-iFlip][0]= max(marker[1-iFlip][0],numberRemaining);
1442          marker[1-iFlip][1]= min(marker[1-iFlip][1],numberPossiblySwapped);
1443          // If we stop now this will be increase in objective (I think)
1444          double increase = (fabs(dualOut_)-totalThru)*theta_;
1445          increase += increaseInObjective;
1446          if (theta_<0.0)
1447            thruThis += fabs(dualOut_); // force using this one
1448          if (increaseInObjective<0.0&&increase<0.0&&lastSequence>=0) {
1449            // back
1450            // We may need to be more careful - we could do by
1451            // switch so we always do fine grained?
1452            bestPivot=0.0;
1453          } else {
1454            // add in
1455            totalThru += thruThis;
1456            increaseInObjective += increaseInThis;
1457          }
1458          if (bestPivot<0.1*bestEverPivot&&
1459              bestEverPivot>1.0e-6&&bestPivot<1.0e-3) {
1460            // back to previous one
1461            sequenceIn_=lastSequence;
1462            // swap regions
1463            iFlip = 1-iFlip;
1464            break;
1465          } else if (sequenceIn_==-1&&upperTheta>largeValue_) {
1466            if (lastPivot>acceptablePivot) {
1467              // back to previous one
1468              sequenceIn_=lastSequence;
1469              // swap regions
1470              iFlip = 1-iFlip;
1471            } else {
1472              // can only get here if all pivots too small
1473            }
1474            break;
1475          } else if (totalThru>=fabs(dualOut_)) {
1476            modifyCosts=true; // fine grain - we can modify costs
1477            break; // no point trying another loop
1478          } else {
1479            lastSequence=sequenceIn_;
1480            if (bestPivot>bestEverPivot)
1481              bestEverPivot=bestPivot;
1482            iFlip = 1 -iFlip;
1483            modifyCosts=true; // fine grain - we can modify costs
1484        }
1485        }
1486        if (iTry==MAXTRY)
1487          iFlip = 1-iFlip; // flip back
1488        break;
1489      } else {
1490        // skip this lot
1491        if (bestPivot>1.0e-3||bestPivot>bestEverPivot) {
1492          bestEverPivot=bestPivot;
1493          lastSequence=bestSequence;
1494        } else {
1495          // keep old swapped
1496          memcpy(array[1-iFlip]+swapped[iFlip],
1497                 array[iFlip]+swapped[iFlip],
1498                 (numberColumns_-swapped[iFlip])*sizeof(double));
1499          memcpy(indices[1-iFlip]+swapped[iFlip],
1500                 indices[iFlip]+swapped[iFlip],
1501                 (numberColumns_-swapped[iFlip])*sizeof(int));
1502          marker[1-iFlip][1] = min(marker[1-iFlip][1],swapped[iFlip]);
1503          swapped[1-iFlip]=swapped[iFlip];
1504        }
1505        increaseInObjective += increaseInThis;
1506        iFlip = 1 - iFlip; // swap regions
1507        tentativeTheta = 2.0*upperTheta;
1508        totalThru += thruThis;
1509      }
1510    }
1511   
1512    // can get here without sequenceIn_ set but with lastSequence
1513    if (sequenceIn_<0&&lastSequence>=0) {
1514      // back to previous one
1515      sequenceIn_=lastSequence;
1516      // swap regions
1517      iFlip = 1-iFlip;
1518    }
1519   
1520    if (sequenceIn_>=0) {
1521      // at this stage sequenceIn_ is just pointer into index array
1522      // flip just so we can use iFlip
1523      iFlip = 1 -iFlip;
1524      spare = array[iFlip];
1525      index = indices[iFlip];
1526      double oldValue;
1527      double alpha = spare[sequenceIn_];
1528      sequenceIn_ = indices[iFlip][sequenceIn_];
1529      oldValue = dj_[sequenceIn_];
1530      theta_ = oldValue/alpha;
1531#if 0
1532      if (numberIterations_>2000)
1533        exit(99);
1534      if (numberIterations_>2000-20)
1535        handler_->setLogLevel(63);
1536      if (numberIterations_>2000-20)
1537        printf("theta %g oldValue %g tol %g %g\n",theta_,oldValue,dualTolerance_,
1538               newTolerance);
1539#endif
1540      if (theta_<MINIMUMTHETA) {
1541        // can't pivot to zero
1542        if (oldValue-MINIMUMTHETA*alpha>=-dualTolerance_) {
1543          theta_=MINIMUMTHETA;
1544        } else if (oldValue-MINIMUMTHETA*alpha>=-newTolerance) {
1545          theta_=MINIMUMTHETA;
1546          thisIncrease=true;
1547        } else {
1548          theta_=(oldValue+newTolerance)/alpha;
1549          assert(theta_>=0.0);
1550          thisIncrease=true;
1551        } 
1552      }
1553      // may need to adjust costs so all dual feasible AND pivoted is exactly 0
1554      if (modifyCosts) {
1555        int i;
1556        double * workRow = rowArray->denseVector();
1557        double * workColumn = columnArray->denseVector();
1558        for (i=numberColumns_-1;i>=swapped[iFlip];i--) {
1559          int iSequence=index[i];
1560          double alpha;
1561          if (iSequence>=numberColumns_)
1562            alpha=workRow[iSequence-numberColumns_];
1563          else
1564            alpha=workColumn[iSequence];
1565          double value = dj_[iSequence]-theta_*alpha;
1566#if 0
1567          if (numberIterations_>2000-20)
1568            printf("%d alpha %g value %g\n",iSequence,alpha,value);
1569#endif
1570           
1571          // can't be free here
1572         
1573          if (alpha < 0.0) {
1574            //at upper bound
1575            if (value>dualTolerance_) {
1576              thisIncrease=true;
1577#define MODIFYCOST 2
1578#if MODIFYCOST
1579              // modify cost to hit new tolerance
1580              double modification = alpha*theta_-dj_[iSequence]
1581                +newTolerance;
1582              //modification = min(modification,dualTolerance_);
1583              //assert (fabs(modification)<1.0e-7);
1584              dj_[iSequence] += modification;
1585              cost_[iSequence] +=  modification;
1586#if 0
1587              if (numberIterations_>2000-20)
1588                printf("%d acost %g mod %g\n",iSequence,cost_[iSequence],
1589                       modification);
1590#endif
1591#endif
1592            }
1593          } else {
1594            // at lower bound
1595            if (-value>dualTolerance_) {
1596              thisIncrease=true;
1597#if MODIFYCOST
1598              // modify cost to hit new tolerance
1599              double modification = alpha*theta_-dj_[iSequence]
1600                -newTolerance;
1601              //modification = max(modification,-dualTolerance_);
1602              //assert (fabs(modification)<1.0e-7);
1603              dj_[iSequence] += modification;
1604              cost_[iSequence] +=  modification;
1605#if 0
1606              if (numberIterations_>2000-20)
1607                printf("%d bcost %g mod %g\n",iSequence,cost_[iSequence],
1608                       modification);
1609#endif
1610#endif
1611            }
1612          }
1613        }
1614      }
1615    }
1616  }
1617
1618  if (sequenceIn_>=0) {
1619    if (sequenceIn_>=numberColumns_) {
1620      //slack
1621      alpha_ = rowArray->denseVector()[sequenceIn_-numberColumns_];
1622    } else {
1623      // column
1624      alpha_ = columnArray->denseVector()[sequenceIn_];
1625    }
1626    lowerIn_ = lower_[sequenceIn_];
1627    upperIn_ = upper_[sequenceIn_];
1628    valueIn_ = solution_[sequenceIn_];
1629    dualIn_ = dj_[sequenceIn_];
1630
1631    if (numberTimesOptimal_) {
1632      // can we adjust cost back closer to original
1633      //*** add coding
1634    }
1635#if MODIFYCOST>1   
1636    // modify cost to hit zero exactly
1637    // so (dualIn_+modification)==theta_*alpha_
1638    double modification = theta_*alpha_-dualIn_;
1639    dualIn_ += modification;
1640    dj_[sequenceIn_]=dualIn_;
1641    cost_[sequenceIn_] += modification;
1642    //assert (fabs(modification)<1.0e-6);
1643#ifdef CLP_DEBUG
1644    if ((handler_->logLevel()&32)&&fabs(modification)>1.0e-15)
1645      printf("exact %d new cost %g, change %g\n",sequenceIn_,
1646             cost_[sequenceIn_],modification);
1647#endif
1648#endif
1649   
1650    if (alpha_<0.0) {
1651      // as if from upper bound
1652      directionIn_=-1;
1653      upperIn_=valueIn_;
1654    } else {
1655      // as if from lower bound
1656      directionIn_=1;
1657      lowerIn_=valueIn_;
1658    }
1659  }
1660  if (thisIncrease) {
1661    newTolerance = dualTolerance_+1.0e-4*dblParam_[ClpDualTolerance];
1662  }
1663
1664  // clear arrays
1665
1666  for (i=0;i<2;i++) {
1667    memset(array[i],0,marker[i][0]*sizeof(double));
1668    memset(array[i]+marker[i][1],0,
1669           (numberColumns_-marker[i][1])*sizeof(double));
1670  }
1671}
1672/* Checks if tentative optimal actually means unbounded
1673   Returns -3 if not, 2 if is unbounded */
1674int 
1675ClpSimplexDual::checkUnbounded(CoinIndexedVector * ray,
1676                                   CoinIndexedVector * spare,
1677                                   double changeCost)
1678{
1679  int status=2; // say unbounded
1680  factorization_->updateColumn(spare,ray);
1681  // get reduced cost
1682  int i;
1683  int number=ray->getNumElements();
1684  int * index = ray->getIndices();
1685  double * array = ray->denseVector();
1686  for (i=0;i<number;i++) {
1687    int iRow=index[i];
1688    int iPivot=pivotVariable_[iRow];
1689    changeCost -= cost(iPivot)*array[iRow];
1690  }
1691  double way;
1692  if (changeCost>0.0) {
1693    //try going down
1694    way=1.0;
1695  } else if (changeCost<0.0) {
1696    //try going up
1697    way=-1.0;
1698  } else {
1699#ifdef CLP_DEBUG
1700    printf("can't decide on up or down\n");
1701#endif
1702    way=0.0;
1703    status=-3;
1704  }
1705  double movement=1.0e10*way; // some largish number
1706  double zeroTolerance = 1.0e-14*dualBound_;
1707  for (i=0;i<number;i++) {
1708    int iRow=index[i];
1709    int iPivot=pivotVariable_[iRow];
1710    double arrayValue = array[iRow];
1711    if (fabs(arrayValue)<zeroTolerance)
1712      arrayValue=0.0;
1713    double newValue=solution(iPivot)+movement*arrayValue;
1714    if (newValue>upper(iPivot)+primalTolerance_||
1715        newValue<lower(iPivot)-primalTolerance_)
1716      status=-3; // not unbounded
1717  }
1718  if (status==2) {
1719    // create ray
1720    delete [] ray_;
1721    ray_ = new double [numberColumns_];
1722    CoinFillN(ray_,numberColumns_,0.0);
1723    for (i=0;i<number;i++) {
1724      int iRow=index[i];
1725      int iPivot=pivotVariable_[iRow];
1726      double arrayValue = array[iRow];
1727      if (iPivot<numberColumns_&&fabs(arrayValue)>=zeroTolerance)
1728        ray_[iPivot] = way* array[iRow];
1729    }
1730  }
1731  ray->clear();
1732  return status;
1733}
1734/* Checks if finished.  Updates status */
1735void 
1736ClpSimplexDual::statusOfProblemInDual(int & lastCleaned,int type,
1737                               ClpSimplexProgress &progress)
1738{
1739  if (type==2) {
1740    // trouble - restore solution
1741    memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
1742    memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
1743           numberRows_*sizeof(double));
1744    memcpy(columnActivityWork_,savedSolution_ ,
1745           numberColumns_*sizeof(double));
1746    forceFactorization_=1; // a bit drastic but ..
1747    changeMade_++; // say something changed
1748  }
1749  int tentativeStatus = problemStatus_;
1750  double changeCost;
1751
1752  if (problemStatus_>-3) {
1753    // factorize
1754    // later on we will need to recover from singularities
1755    // also we could skip if first time
1756    // save dual weights
1757    dualRowPivot_->saveWeights(this,1);
1758    if (type) {
1759      // is factorization okay?
1760      if (internalFactorize(1)) {
1761        // no - restore previous basis
1762        assert (type==1);
1763        changeMade_++; // say something changed
1764        memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
1765        memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
1766               numberRows_*sizeof(double));
1767        memcpy(columnActivityWork_,savedSolution_ ,
1768               numberColumns_*sizeof(double));
1769        // get correct bounds on all variables
1770        double dummyChangeCost=0.0;
1771        changeBounds(true,rowArray_[2],dummyChangeCost);
1772        // throw away change
1773        rowArray_[2]->clear();
1774        forceFactorization_=1; // a bit drastic but ..
1775        type = 2;
1776        assert (internalFactorize(1)==0);
1777      }
1778    }
1779    problemStatus_=-3;
1780  }
1781  // at this stage status is -3 or -4 if looks infeasible
1782  // get primal and dual solutions
1783  gutsOfSolution(rowActivityWork_,columnActivityWork_);
1784  // Check if looping
1785  int loop = progress.looping();
1786  bool situationChanged=false;
1787  if (loop>=0) {
1788    problemStatus_ = loop; //exit if in loop
1789    return;
1790  } else if (loop<-1) {
1791    // something may have changed
1792    gutsOfSolution(rowActivityWork_,columnActivityWork_);
1793    situationChanged=true;
1794  }
1795  progressFlag_ = 0; //reset progress flag
1796#ifdef CLP_DEBUG
1797  if (!rowScale_&&(handler_->logLevel()&32)) {
1798    double * objectiveSimplex
1799      = copyOfArray(objective_,numberColumns_,0.0);
1800    double * rowObjectiveSimplex
1801      = copyOfArray(rowObjective_,numberRows_,0.0);
1802    int i;
1803    double largest;
1804    largest=0.0;
1805    for (i=0;i<numberRows_;i++) {
1806      rowObjectiveSimplex[i] *= optimizationDirection_;
1807      double difference = fabs(rowObjectiveWork_[i]-rowObjectiveSimplex[i]);
1808      if (difference>largest)
1809        largest=difference;
1810    }
1811    for (i=0;i<numberColumns_;i++) {
1812      objectiveSimplex[i] *= optimizationDirection_;
1813      double difference = fabs(objectiveWork_[i]-objectiveSimplex[i]);
1814      if (difference>largest)
1815        largest=difference;
1816    }
1817    if ((handler_->logLevel()&16))
1818      printf("difference in obj %g\n",largest);
1819    delete [] objectiveSimplex;
1820    delete [] rowObjectiveSimplex;
1821  }
1822#endif
1823  handler_->message(CLP_SIMPLEX_STATUS,messages_)
1824    <<numberIterations_<<objectiveValue();
1825  handler_->printing(sumPrimalInfeasibilities_>0.0)
1826    <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
1827  handler_->printing(sumDualInfeasibilities_>0.0)
1828    <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
1829  handler_->printing(numberDualInfeasibilitiesWithoutFree_
1830                     <numberDualInfeasibilities_)
1831                       <<numberDualInfeasibilities_-
1832    numberDualInfeasibilitiesWithoutFree_;
1833  handler_->message()<<CoinMessageEol;
1834
1835  while (problemStatus_<=-3) {
1836    bool cleanDuals=false;
1837    int numberChangedBounds=0;
1838    int doOriginalTolerance=0;
1839    if ( lastCleaned==numberIterations_)
1840      doOriginalTolerance=1;
1841    // check optimal
1842    // give code benefit of doubt
1843    if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
1844        sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
1845      // say optimal (with these bounds etc)
1846      numberDualInfeasibilities_ = 0;
1847      sumDualInfeasibilities_ = 0.0;
1848      numberPrimalInfeasibilities_ = 0;
1849      sumPrimalInfeasibilities_ = 0.0;
1850    }
1851    if (dualFeasible()||problemStatus_==-4) {
1852      if (primalFeasible()) {
1853        // may be optimal - or may be bounds are wrong
1854        handler_->message(CLP_DUAL_BOUNDS,messages_)
1855          <<dualBound_
1856          <<CoinMessageEol;
1857        // save solution in case unbounded
1858        CoinDisjointCopyN(columnActivityWork_,numberColumns_,
1859                          columnArray_[0]->denseVector());
1860        CoinDisjointCopyN(rowActivityWork_,numberRows_,
1861                          rowArray_[2]->denseVector());
1862        numberChangedBounds=changeBounds(false,rowArray_[0],changeCost);
1863        if (numberChangedBounds<=0) {
1864          //looks optimal - do we need to reset tolerance
1865          if (lastCleaned<numberIterations_&&numberTimesOptimal_<4) {
1866            doOriginalTolerance=2;
1867            numberTimesOptimal_++;
1868            changeMade_++; // say something changed
1869            if (numberTimesOptimal_==1) {
1870              dualTolerance_ = min(dualTolerance_,1.0e-8);
1871              // better to have small tolerance even if slower
1872              factorization_->zeroTolerance(1.0e-15);
1873            }
1874            cleanDuals=true;
1875          } else {
1876            problemStatus_=0; // optimal
1877            if (lastCleaned<numberIterations_) {
1878              handler_->message(CLP_SIMPLEX_GIVINGUP,messages_)
1879                <<CoinMessageEol;
1880            }
1881          }
1882        } else {
1883          cleanDuals=true;
1884          if (doOriginalTolerance==1) {
1885            // check unbounded
1886            problemStatus_ = checkUnbounded(rowArray_[0],rowArray_[1],
1887                                            changeCost);
1888            if (problemStatus_==2) {
1889              // it is unbounded - restore solution
1890              // but first add in changes to non-basic
1891              int iColumn;
1892              double * original = columnArray_[0]->denseVector();
1893              for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1894                if(getColumnStatus(iColumn)!= ClpSimplex::basic)
1895                  ray_[iColumn] += 
1896                    columnActivityWork_[iColumn]-original[iColumn];
1897                columnActivityWork_[iColumn] = original[iColumn];
1898              }
1899              CoinDisjointCopyN(rowArray_[2]->denseVector(),numberRows_,
1900                                rowActivityWork_);
1901            }
1902          } else {
1903            doOriginalTolerance=2;
1904            rowArray_[0]->clear();
1905          }
1906        }
1907        CoinFillN(columnArray_[0]->denseVector(),numberColumns_,0.0);
1908        CoinFillN(rowArray_[2]->denseVector(),numberRows_,0.0);
1909      } 
1910      if (problemStatus_==-4) {
1911        // may be infeasible - or may be bounds are wrong
1912        handler_->message(CLP_DUAL_CHECKB,messages_)
1913          <<dualBound_
1914          <<CoinMessageEol;
1915        numberChangedBounds=changeBounds(false,NULL,changeCost);
1916        if (numberChangedBounds<=0||dualBound_>1.0e20||
1917            (largestPrimalError_>1.0&&dualBound_>1.0e17)) {
1918          problemStatus_=1; // infeasible
1919        } else {
1920          problemStatus_=-1; //iterate
1921          cleanDuals=true;
1922          doOriginalTolerance=2;
1923          // and delete ray which has been created
1924          delete [] ray_;
1925          ray_ = NULL;
1926        }
1927      }
1928    } else {
1929      cleanDuals=true;
1930    }
1931    if (problemStatus_<0) {
1932      if (doOriginalTolerance==2) {
1933        // put back original tolerance
1934        lastCleaned=numberIterations_;
1935        handler_->message(CLP_DUAL_ORIGINAL,messages_)
1936          <<CoinMessageEol;
1937
1938        perturbation_=102; // stop any perturbations
1939        createRim(4);
1940        // make sure duals are current
1941        computeDuals();
1942        // put back bounds as they were if was optimal
1943        if (doOriginalTolerance==2) {
1944          changeMade_++; // say something changed
1945          changeBounds(true,NULL,changeCost);
1946          cleanDuals=true;
1947        }
1948      }
1949      if (cleanDuals) {
1950        // make sure dual feasible
1951        // look at all rows and columns
1952        rowArray_[0]->clear();
1953        CoinIotaN(rowArray_[0]->getIndices(),numberRows_,0);
1954        rowArray_[0]->setNumElements(numberRows_);
1955        columnArray_[0]->clear();
1956        CoinIotaN(columnArray_[0]->getIndices(),numberColumns_,0);
1957        columnArray_[0]->setNumElements(numberColumns_);
1958        double objectiveChange=0.0;
1959        updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[1],
1960                          0.0,objectiveChange);
1961        // for now - recompute all
1962        gutsOfSolution(rowActivityWork_,columnActivityWork_);
1963        assert(numberDualInfeasibilitiesWithoutFree_==0);
1964        if (numberDualInfeasibilities_) {
1965          // bad free variables
1966          if (primalFeasible()) {
1967            std::cerr<<"Free variable problem?"<<std::endl;
1968            abort(); // what now
1969          }
1970          problemStatus_=-1; // carry on as normal
1971        }
1972      } else {
1973        // iterate
1974        problemStatus_=-1;
1975      }
1976    }
1977  }
1978  if (type==0||type==1) {
1979    if (!type) {
1980      // create save arrays
1981      delete [] saveStatus_;
1982      delete [] savedSolution_;
1983      saveStatus_ = new unsigned char [numberRows_+numberColumns_];
1984      savedSolution_ = new double [numberRows_+numberColumns_];
1985    }
1986    // save arrays
1987    memcpy(saveStatus_,status_,(numberColumns_+numberRows_)*sizeof(char));
1988    memcpy(savedSolution_+numberColumns_ ,rowActivityWork_,
1989           numberRows_*sizeof(double));
1990    memcpy(savedSolution_ ,columnActivityWork_,numberColumns_*sizeof(double));
1991  }
1992
1993  // restore weights (if saved) - also recompute infeasibility list
1994  if (tentativeStatus>-3) 
1995    dualRowPivot_->saveWeights(this,(type <2) ? 2 : 4);
1996  else
1997    dualRowPivot_->saveWeights(this,3);
1998  // unflag all variables (we may want to wait a bit?)
1999  int iRow;
2000  for (iRow=0;iRow<numberRows_;iRow++) {
2001    int iPivot=pivotVariable_[iRow];
2002    clearFlagged(iPivot);
2003  }
2004  if (problemStatus_<0&&!changeMade_) {
2005    problemStatus_=4; // unknown
2006  }
2007
2008}
2009/* While updateDualsInDual sees what effect is of flip
2010   this does actual flipping.
2011   If change >0.0 then value in array >0.0 => from lower to upper
2012*/
2013void 
2014ClpSimplexDual::flipBounds(CoinIndexedVector * rowArray,
2015                  CoinIndexedVector * columnArray,
2016                  double change)
2017{
2018  double * work;
2019  int number;
2020  int * which;
2021 
2022  int iSection;
2023
2024  for (iSection=0;iSection<2;iSection++) {
2025    int i;
2026    double * solution = solutionRegion(iSection);
2027    double * lower = lowerRegion(iSection);
2028    double * upper = upperRegion(iSection);
2029    int addSequence;
2030    if (!iSection) {
2031      work = rowArray->denseVector();
2032      number = rowArray->getNumElements();
2033      which = rowArray->getIndices();
2034      addSequence = numberColumns_;
2035    } else {
2036      work = columnArray->denseVector();
2037      number = columnArray->getNumElements();
2038      which = columnArray->getIndices();
2039      addSequence = 0;
2040    }
2041   
2042    for (i=0;i<number;i++) {
2043      int iSequence = which[i];
2044#ifndef NDEBUG
2045      double value = work[iSequence]*change;
2046#endif
2047      work[iSequence]=0.0;
2048      Status status = getStatus(iSequence+addSequence);
2049
2050      switch(status) {
2051
2052      case ClpSimplex::basic:
2053      case ClpSimplex::isFree:
2054      case ClpSimplex::superBasic:
2055        break;
2056      case ClpSimplex::atUpperBound:
2057        assert (value<=0.0);
2058        // to lower bound
2059        setStatus(iSequence+addSequence,ClpSimplex::atLowerBound);
2060        solution[iSequence] = lower[iSequence];
2061        break;
2062      case ClpSimplex::atLowerBound:
2063        assert (value>=0.0);
2064        // to upper bound
2065        setStatus(iSequence+addSequence,ClpSimplex::atUpperBound);
2066        solution[iSequence] = upper[iSequence];
2067        break;
2068      }
2069    }
2070  }
2071  rowArray->setNumElements(0);
2072  columnArray->setNumElements(0);
2073}
2074// Restores bound to original bound
2075void 
2076ClpSimplexDual::originalBound( int iSequence)
2077{
2078  if (iSequence>=numberColumns_) {
2079    // rows
2080    int iRow = iSequence-numberColumns_;
2081    rowLowerWork_[iRow]=rowLower_[iRow];
2082    rowUpperWork_[iRow]=rowUpper_[iRow];
2083    if (rowScale_) {
2084      if (rowLowerWork_[iRow]>-1.0e50)
2085        rowLowerWork_[iRow] *= rowScale_[iRow];
2086      if (rowUpperWork_[iRow]<1.0e50)
2087        rowUpperWork_[iRow] *= rowScale_[iRow];
2088    }
2089  } else {
2090    // columns
2091    columnLowerWork_[iSequence]=columnLower_[iSequence];
2092    columnUpperWork_[iSequence]=columnUpper_[iSequence];
2093    if (rowScale_) {
2094      double multiplier = 1.0/columnScale_[iSequence];
2095      if (columnLowerWork_[iSequence]>-1.0e50)
2096        columnLowerWork_[iSequence] *= multiplier;
2097      if (columnUpperWork_[iSequence]<1.0e50)
2098        columnUpperWork_[iSequence] *= multiplier;
2099    }
2100  }
2101  setFakeBound(iSequence,ClpSimplex::noFake);
2102}
2103/* As changeBounds but just changes new bounds for a single variable.
2104   Returns true if change */
2105bool 
2106ClpSimplexDual::changeBound( int iSequence)
2107{
2108  // old values
2109  double oldLower=lower_[iSequence];
2110  double oldUpper=upper_[iSequence];
2111  double value=solution_[iSequence];
2112  bool modified=false;
2113  originalBound(iSequence);
2114  // original values
2115  double lowerValue=lower_[iSequence];
2116  double upperValue=upper_[iSequence];
2117  // back to altered values
2118  lower_[iSequence] = oldLower;
2119  upper_[iSequence] = oldUpper;
2120  if (value==oldLower) {
2121    if (upperValue > oldLower + dualBound_) {
2122      upper_[iSequence]=oldLower+dualBound_;
2123      setFakeBound(iSequence,ClpSimplex::upperFake);
2124      modified=true;
2125    }
2126  } else if (value==oldUpper) {
2127    if (lowerValue < oldUpper - dualBound_) {
2128      lower_[iSequence]=oldUpper-dualBound_;
2129      setFakeBound(iSequence,ClpSimplex::lowerFake);
2130      modified=true;
2131    }
2132  } else {
2133    assert(value==oldLower||value==oldUpper);
2134  }
2135  return modified;
2136}
2137// Perturbs problem
2138void 
2139ClpSimplexDual::perturb()
2140{
2141  if (perturbation_>100)
2142    return; //perturbed already
2143  int iRow,iColumn;
2144  // dual perturbation
2145  double perturbation=1.0e-20;
2146  // maximum fraction of cost to perturb
2147  double maximumFraction = 1.0e-4;
2148  if (perturbation_==100) {
2149    perturbation = 1.0e-4;
2150    for (iRow=0;iRow<numberRows_;iRow++) {
2151      double value = fabs(rowActivityWork_[iRow]*rowObjectiveWork_[iRow]);
2152      perturbation = max(perturbation,value);
2153    }
2154    for (iColumn=0;iColumn<numberColumns_;iColumn++) { 
2155      double value = 
2156        fabs(columnActivityWork_[iColumn]*objectiveWork_[iColumn]);
2157      perturbation = max(perturbation,value);
2158    }
2159    perturbation *= 1.0e-8;
2160  } else if (perturbation_<100) {
2161    perturbation = pow(10.0,perturbation_);
2162    // user is in charge
2163    maximumFraction = 1.0e100;
2164  }
2165  // modify costs
2166  handler_->message(CLP_SIMPLEX_PERTURB,messages_)
2167    <<perturbation
2168    <<CoinMessageEol;
2169  for (iRow=0;iRow<numberRows_;iRow++) {
2170    double value = perturbation;
2171    double currentValue = rowObjectiveWork_[iRow];
2172    value = min(value,maximumFraction*fabs(currentValue)+1.0e-6);
2173    if (rowLowerWork_[iRow]>-largeValue_) {
2174      if (fabs(rowLowerWork_[iRow])<fabs(rowUpperWork_[iRow])) 
2175        value *= drand48();
2176      else
2177        value *= -drand48();
2178    } else if (rowUpperWork_[iRow]<largeValue_) {
2179      value *= -drand48();
2180    } else {
2181      value=0.0;
2182    }
2183    rowObjectiveWork_[iRow] += value;
2184  }
2185  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
2186    double value = perturbation;
2187    double currentValue = objectiveWork_[iColumn];
2188    value = min(value,maximumFraction*fabs(currentValue)+1.0e-6);
2189    if (columnLowerWork_[iColumn]>-largeValue_) {
2190      if (fabs(columnLowerWork_[iColumn])<
2191          fabs(columnUpperWork_[iColumn])) 
2192        value *= drand48();
2193      else
2194        value *= -drand48();
2195    } else if (columnUpperWork_[iColumn]<largeValue_) {
2196      value *= -drand48();
2197    } else {
2198      value=0.0;
2199    }
2200    objectiveWork_[iColumn] += value;
2201  }
2202  // say perturbed
2203  perturbation_=102;
2204
2205}
2206/* For strong branching.  On input lower and upper are new bounds
2207   while on output they are change in objective function values
2208   (>1.0e50 infeasible).
2209   Return code is 0 if nothing interesting, -1 if infeasible both
2210   ways and +1 if infeasible one way (check values to see which one(s))
2211*/
2212int ClpSimplexDual::strongBranching(int numberVariables,const int * variables,
2213                                    double * newLower, double * newUpper,
2214                                    bool stopOnFirstInfeasible,
2215                                    bool alwaysFinish)
2216{
2217  int i;
2218  int returnCode=0;
2219  double saveObjectiveValue = objectiveValue_;
2220#if 1
2221  algorithm_ = -1;
2222
2223  //scaling(false);
2224
2225  // put in standard form (and make row copy)
2226  // create modifiable copies of model rim and do optional scaling
2227  createRim(7+8+16,true);
2228
2229  // change newLower and newUpper if scaled
2230
2231  // Do initial factorization
2232  // and set certain stuff
2233  // We can either set increasing rows so ...IsBasic gives pivot row
2234  // or we can just increment iBasic one by one
2235  // for now let ...iBasic give pivot row
2236  factorization_->increasingRows(2);
2237  // row activities have negative sign
2238  factorization_->slackValue(-1.0);
2239  factorization_->zeroTolerance(1.0e-13);
2240  // save if sparse factorization wanted
2241  int saveSparse = factorization_->sparseThreshold();
2242
2243  int factorizationStatus = internalFactorize(0);
2244  if (factorizationStatus<0)
2245    return 1; // some error
2246  else if (factorizationStatus)
2247    handler_->message(CLP_SINGULARITIES,messages_)
2248      <<factorizationStatus
2249      <<CoinMessageEol;
2250  if (saveSparse) {
2251    // use default at present
2252    factorization_->sparseThreshold(0);
2253    factorization_->goSparse();
2254  }
2255 
2256  // save stuff
2257  ClpFactorization saveFactorization(*factorization_);
2258  // save basis and solution
2259  double * saveSolution = new double[numberRows_+numberColumns_];
2260  memcpy(saveSolution,solution_,
2261         (numberRows_+numberColumns_)*sizeof(double));
2262  unsigned char * saveStatus = 
2263    new unsigned char [numberRows_+numberColumns_];
2264  memcpy(saveStatus,status_,(numberColumns_+numberRows_)*sizeof(char));
2265  // save bounds as createRim makes clean copies
2266  double * saveLower = new double[numberRows_+numberColumns_];
2267  memcpy(saveLower,lower_,
2268         (numberRows_+numberColumns_)*sizeof(double));
2269  double * saveUpper = new double[numberRows_+numberColumns_];
2270  memcpy(saveUpper,upper_,
2271         (numberRows_+numberColumns_)*sizeof(double));
2272  double * saveObjective = new double[numberRows_+numberColumns_];
2273  memcpy(saveObjective,cost_,
2274         (numberRows_+numberColumns_)*sizeof(double));
2275  int * savePivot = new int [numberRows_];
2276  memcpy(savePivot, pivotVariable_, numberRows_*sizeof(int));
2277  // need to save/restore weights.
2278
2279  for (i=0;i<numberVariables;i++) {
2280    int iColumn = variables[i];
2281    double objectiveChange;
2282    double saveBound;
2283   
2284    // try down
2285   
2286    saveBound = columnUpper_[iColumn];
2287    // external view - in case really getting optimal
2288    columnUpper_[iColumn] = newUpper[i];
2289    if (scalingFlag_<=0) 
2290      upper_[iColumn] = newUpper[i];
2291    else 
2292      upper_[iColumn] = newUpper[i]/columnScale_[iColumn]; // scale
2293    // Start of fast iterations
2294    int status = fastDual(alwaysFinish);
2295
2296    // restore
2297    memcpy(solution_,saveSolution,
2298           (numberRows_+numberColumns_)*sizeof(double));
2299    memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char));
2300    memcpy(lower_,saveLower,
2301           (numberRows_+numberColumns_)*sizeof(double));
2302    memcpy(upper_,saveUpper,
2303           (numberRows_+numberColumns_)*sizeof(double));
2304    memcpy(cost_,saveObjective,
2305         (numberRows_+numberColumns_)*sizeof(double));
2306    columnUpper_[iColumn] = saveBound;
2307    memcpy(pivotVariable_, savePivot, numberRows_*sizeof(int));
2308    delete factorization_;
2309    factorization_ = new ClpFactorization(saveFactorization);
2310
2311    if (status||(problemStatus_==0&&!isDualObjectiveLimitReached())) {
2312      objectiveChange = objectiveValue_-saveObjectiveValue;
2313    } else {
2314      objectiveChange = 1.0e100;
2315    }
2316    newUpper[i]=objectiveChange;
2317#ifdef CLP_DEBUG
2318    printf("down on %d costs %g\n",iColumn,objectiveChange);
2319#endif
2320         
2321    // try up
2322   
2323    saveBound = columnLower_[iColumn];
2324    // external view - in case really getting optimal
2325    columnLower_[iColumn] = newLower[i];
2326    if (scalingFlag_<=0) 
2327      lower_[iColumn] = newLower[i];
2328    else 
2329      lower_[iColumn] = newLower[i]/columnScale_[iColumn]; // scale
2330    // Start of fast iterations
2331    status = fastDual(alwaysFinish);
2332
2333    // restore
2334    memcpy(solution_,saveSolution,
2335           (numberRows_+numberColumns_)*sizeof(double));
2336    memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char));
2337    memcpy(lower_,saveLower,
2338           (numberRows_+numberColumns_)*sizeof(double));
2339    memcpy(upper_,saveUpper,
2340           (numberRows_+numberColumns_)*sizeof(double));
2341    memcpy(cost_,saveObjective,
2342         (numberRows_+numberColumns_)*sizeof(double));
2343    columnLower_[iColumn] = saveBound;
2344    memcpy(pivotVariable_, savePivot, numberRows_*sizeof(int));
2345    delete factorization_;
2346    factorization_ = new ClpFactorization(saveFactorization);
2347
2348    if (status||(problemStatus_==0&&!isDualObjectiveLimitReached())) {
2349      objectiveChange = objectiveValue_-saveObjectiveValue;
2350    } else {
2351      objectiveChange = 1.0e100;
2352    }
2353    newLower[i]=objectiveChange;
2354#ifdef CLP_DEBUG
2355    printf("up on %d costs %g\n",iColumn,objectiveChange);
2356#endif
2357         
2358    /* Possibilities are:
2359       Both sides feasible - store
2360       Neither side feasible - set objective high and exit
2361       One side feasible - change bounds and resolve
2362    */
2363    if (newUpper[i]<1.0e100) {
2364      if(newLower[i]<1.0e100) {
2365        // feasible - no action
2366      } else {
2367        // up feasible, down infeasible
2368        returnCode=1;
2369        if (stopOnFirstInfeasible)
2370          break;
2371      }
2372    } else {
2373      if(newLower[i]<1.0e100) {
2374        // down feasible, up infeasible
2375        returnCode=1;
2376        if (stopOnFirstInfeasible)
2377          break;
2378      } else {
2379        // neither side feasible
2380        returnCode=-1;
2381        break;
2382      }
2383    }
2384  }
2385  delete [] saveSolution;
2386  delete [] saveLower;
2387  delete [] saveUpper;
2388  delete [] saveObjective;
2389  delete [] saveStatus;
2390  delete [] savePivot;
2391
2392  // at present we are leaving factorization around
2393  // maybe we should empty it
2394  deleteRim();
2395  factorization_->sparseThreshold(saveSparse);
2396#else
2397  // save basis and solution
2398  double * rowSolution = new double[numberRows_];
2399  memcpy(rowSolution,rowActivity_,numberRows_*sizeof(double));
2400  double * columnSolution = new double[numberColumns_];
2401  memcpy(columnSolution,columnActivity_,numberColumns_*sizeof(double));
2402  unsigned char * saveStatus = 
2403    new unsigned char [numberRows_+numberColumns_];
2404  if (!status_) {
2405    // odd but anyway setup all slack basis
2406    status_ = new unsigned char [numberColumns_+numberRows_];
2407    memset(status_,0,(numberColumns_+numberRows_)*sizeof(char));
2408  }
2409  memcpy(saveStatus,status_,(numberColumns_+numberRows_)*sizeof(char));
2410  for (i=0;i<numberVariables;i++) {
2411    int iColumn = variables[i];
2412    double objectiveChange;
2413   
2414    // try down
2415   
2416    double saveUpper = columnUpper_[iColumn];
2417    columnUpper_[iColumn] = newUpper[i];
2418    dual();
2419
2420    // restore
2421    columnUpper_[iColumn] = saveUpper;
2422    memcpy(rowActivity_,rowSolution,numberRows_*sizeof(double));
2423    memcpy(columnActivity_,columnSolution,numberColumns_*sizeof(double));
2424    memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char));
2425
2426    if (problemStatus_==0&&!isDualObjectiveLimitReached()) {
2427      objectiveChange = objectiveValue_-saveObjectiveValue;
2428    } else {
2429      objectiveChange = 1.0e100;
2430    }
2431    newUpper[i]=objectiveChange;
2432#ifdef CLP_DEBUG
2433    printf("down on %d costs %g\n",iColumn,objectiveChange);
2434#endif
2435         
2436    // try up
2437   
2438    double saveLower = columnLower_[iColumn];
2439    columnLower_[iColumn] = newLower[i];
2440    dual();
2441
2442    // restore
2443    columnLower_[iColumn] = saveLower;
2444    memcpy(rowActivity_,rowSolution,numberRows_*sizeof(double));
2445    memcpy(columnActivity_,columnSolution,numberColumns_*sizeof(double));
2446    memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char));
2447
2448    if (problemStatus_==0&&!isDualObjectiveLimitReached()) {
2449      objectiveChange = objectiveValue_-saveObjectiveValue;
2450    } else {
2451      objectiveChange = 1.0e100;
2452    }
2453    newLower[i]=objectiveChange;
2454#ifdef CLP_DEBUG
2455    printf("up on %d costs %g\n",iColumn,objectiveChange);
2456#endif
2457         
2458    /* Possibilities are:
2459       Both sides feasible - store
2460       Neither side feasible - set objective high and exit
2461       One side feasible - change bounds and resolve
2462    */
2463    if (newUpper[i]<1.0e100) {
2464      if(newLower[i]<1.0e100) {
2465        // feasible - no action
2466      } else {
2467        // up feasible, down infeasible
2468        returnCode=1;
2469        if (stopOnFirstInfeasible)
2470          break;
2471      }
2472    } else {
2473      if(newLower[i]<1.0e100) {
2474        // down feasible, up infeasible
2475        returnCode=1;
2476        if (stopOnFirstInfeasible)
2477          break;
2478      } else {
2479        // neither side feasible
2480        returnCode=-1;
2481        break;
2482      }
2483    }
2484  }
2485  delete [] rowSolution;
2486  delete [] columnSolution;
2487  delete [] saveStatus;
2488#endif
2489  objectiveValue_ = saveObjectiveValue;
2490  return returnCode;
2491}
2492// treat no pivot as finished (unless interesting)
2493int ClpSimplexDual::fastDual(bool alwaysFinish)
2494{
2495  algorithm_ = -1;
2496  dualTolerance_=dblParam_[ClpDualTolerance];
2497  primalTolerance_=dblParam_[ClpPrimalTolerance];
2498
2499  // save dual bound
2500  double saveDualBound = dualBound_;
2501
2502  double objectiveChange;
2503  // for dual we will change bounds using dualBound_
2504  // for this we need clean basis so it is after factorize
2505  gutsOfSolution(rowActivityWork_,columnActivityWork_);
2506  numberFake_ =0; // Number of variables at fake bounds
2507  changeBounds(true,NULL,objectiveChange);
2508
2509  problemStatus_ = -1;
2510  numberIterations_=0;
2511
2512  int lastCleaned=0; // last time objective or bounds cleaned up
2513
2514  // number of times we have declared optimality
2515  numberTimesOptimal_=0;
2516
2517  // Progress indicator
2518  ClpSimplexProgress progress(this);
2519
2520  // This says whether to restore things etc
2521  int factorType=0;
2522  /*
2523    Status of problem:
2524    0 - optimal
2525    1 - infeasible
2526    2 - unbounded
2527    -1 - iterating
2528    -2 - factorization wanted
2529    -3 - redo checking without factorization
2530    -4 - looks infeasible
2531
2532    BUT also from whileIterating return code is:
2533
2534   -1 iterations etc
2535   -2 inaccuracy
2536   -3 slight inaccuracy (and done iterations)
2537   +0 looks optimal (might be unbounded - but we will investigate)
2538   +1 looks infeasible
2539   +3 max iterations
2540
2541  */
2542
2543  int returnCode = 0;
2544
2545  while (problemStatus_<0) {
2546    int iRow,iColumn;
2547    // clear
2548    for (iRow=0;iRow<4;iRow++) {
2549      rowArray_[iRow]->clear();
2550    }   
2551   
2552    for (iColumn=0;iColumn<2;iColumn++) {
2553      columnArray_[iColumn]->clear();
2554    }   
2555
2556    // give matrix (and model costs and bounds a chance to be
2557    // refreshed (normally null)
2558    matrix_->refresh(this);
2559    // may factorize, checks if problem finished
2560    // should be able to speed this up on first time
2561    statusOfProblemInDual(lastCleaned,factorType,progress);
2562
2563    // Say good factorization
2564    factorType=1;
2565
2566    // Do iterations
2567    if (problemStatus_<0) {
2568      returnCode = whileIterating();
2569      if (!alwaysFinish&&returnCode<1) {
2570        double limit = 0.0;
2571        getDblParam(ClpDualObjectiveLimit, limit);
2572        if(objectiveValue()*optimizationDirection_<limit|| 
2573           numberAtFakeBound()) {
2574          // can't say anything interesting - might as well return
2575#ifdef CLP_DEBUG
2576          printf("returning from fastDual after %d iterations with code %d\n",
2577                 numberIterations_,returnCode);
2578#endif
2579          returnCode=1;
2580          break;
2581        }
2582      }
2583      returnCode=0;
2584    }
2585  }
2586
2587  assert(!numberFake_||returnCode||problemStatus_); // all bounds should be okay
2588  dualBound_ = saveDualBound;
2589  return returnCode;
2590}
2591/* Checks number of variables at fake bounds.  This is used by fastDual
2592   so can exit gracefully before end */
2593int 
2594ClpSimplexDual::numberAtFakeBound()
2595{
2596  int iSequence;
2597  int numberFake=0;
2598 
2599  for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
2600    FakeBound bound = getFakeBound(iSequence);
2601    switch(getStatus(iSequence)) {
2602
2603    case ClpSimplex::basic:
2604      break;
2605    case ClpSimplex::isFree:
2606    case ClpSimplex::superBasic:
2607      break;
2608    case ClpSimplex::atUpperBound:
2609      if (bound==ClpSimplex::upperFake||bound==ClpSimplex::bothFake)
2610        numberFake++;
2611      break;
2612    case ClpSimplex::atLowerBound:
2613      if (bound==ClpSimplex::lowerFake||bound==ClpSimplex::bothFake)
2614        numberFake++;
2615      break;
2616    }
2617  }
2618  numberFake_ = numberFake;
2619  return numberFake;
2620}
Note: See TracBrowser for help on using the repository browser.