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

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

For presolve

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