source: trunk/ClpSimplexDual.cpp @ 225

Last change on this file since 225 was 225, checked in by forrest, 16 years ago

This should break everything

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 120.0 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
90#include "CoinPragma.hpp"
91
92#include <math.h>
93
94#include "CoinHelperFunctions.hpp"
95#include "ClpSimplexDual.hpp"
96#include "ClpFactorization.hpp"
97#include "CoinPackedMatrix.hpp"
98#include "CoinIndexedVector.hpp"
99#include "CoinWarmStartBasis.hpp"
100#include "ClpDualRowDantzig.hpp"
101#include "ClpPlusMinusOneMatrix.hpp"
102#include "ClpMessage.hpp"
103#include <cfloat>
104#include <cassert>
105#include <string>
106#include <stdio.h>
107#include <iostream>
108//#define CLP_DEBUG 1
109#define CHECK_DJ
110// dual
111int ClpSimplexDual::dual (int ifValuesPass )
112{
113
114  /* *** Method
115     This is a vanilla version of dual simplex.
116
117     It tries to be a single phase approach with a weight of 1.0 being
118     given to getting optimal and a weight of dualBound_ being
119     given to getting dual feasible.  In this version I have used the
120     idea that this weight can be thought of as a fake bound.  If the
121     distance between the lower and upper bounds on a variable is less
122     than the feasibility weight then we are always better off flipping
123     to other bound to make dual feasible.  If the distance is greater
124     then we make up a fake bound dualBound_ away from one bound.
125     If we end up optimal or primal infeasible, we check to see if
126     bounds okay.  If so we have finished, if not we increase dualBound_
127     and continue (after checking if unbounded). I am undecided about
128     free variables - there is coding but I am not sure about it.  At
129     present I put them in basis anyway.
130
131     The code is designed to take advantage of sparsity so arrays are
132     seldom zeroed out from scratch or gone over in their entirety.
133     The only exception is a full scan to find outgoing variable.  This
134     will be changed to keep an updated list of infeasibilities (or squares
135     if steepest edge).  Also on easy problems we don't need full scan - just
136     pick first reasonable.
137
138     One problem is how to tackle degeneracy and accuracy.  At present
139     I am using the modification of costs which I put in OSL and which was
140     extended by Gill et al.  I am still not sure of the exact details.
141
142     The flow of dual is three while loops as follows:
143
144     while (not finished) {
145
146       while (not clean solution) {
147
148          Factorize and/or clean up solution by flipping variables so
149          dual feasible.  If looks finished check fake dual bounds.
150          Repeat until status is iterating (-1) or finished (0,1,2)
151
152       }
153
154       while (status==-1) {
155
156         Iterate until no pivot in or out or time to re-factorize.
157
158         Flow is:
159
160         choose pivot row (outgoing variable).  if none then
161         we are primal feasible so looks as if done but we need to
162         break and check bounds etc.
163
164         Get pivot row in tableau
165
166         Choose incoming column.  If we don't find one then we look
167         primal infeasible so break and check bounds etc.  (Also the
168         pivot tolerance is larger after any iterations so that may be
169         reason)
170
171         If we do find incoming column, we may have to adjust costs to
172         keep going forwards (anti-degeneracy).  Check pivot will be stable
173         and if unstable throw away iteration (we will need to implement
174         flagging of basic variables sometime) and break to re-factorize.
175         If minor error re-factorize after iteration.
176
177         Update everything (this may involve flipping variables to stay
178         dual feasible.
179
180       }
181
182     }
183
184     At present we never check we are going forwards.  I overdid that in
185     OSL so will try and make a last resort.
186
187     Needs partial scan pivot out option.
188     Needs dantzig, uninitialized and full steepest edge options (can still
189     use partial scan)
190
191     May need other anti-degeneracy measures, especially if we try and use
192     loose tolerances as a way to solve in fewer iterations.
193
194     I like idea of dynamic scaling.  This gives opportunity to decouple
195     different implications of scaling for accuracy, iteration count and
196     feasibility tolerance.
197
198  */
199  algorithm_ = -1;
200
201  // save data
202  ClpDataSave data = saveData();
203
204  // Save so can see if doing after primal
205  int initialStatus=problemStatus_;
206
207  // If values pass then save given duals round check solution
208  double * saveDuals = NULL;
209  if (ifValuesPass) {
210    saveDuals = new double [numberRows_+numberColumns_];
211    memcpy(saveDuals,dual_,numberRows_*sizeof(double));
212  }
213 
214  // sanity check
215  // initialize - no values pass and algorithm_ is -1
216  // put in standard form (and make row copy)
217  // create modifiable copies of model rim and do optional scaling
218  // If problem looks okay
219  // Do initial factorization
220  // If user asked for perturbation - do it
221  if (!startup(0)) {
222    // looks okay
223   
224    // If values pass then scale pi
225    if (ifValuesPass) {
226      if (problemStatus_&&perturbation_<100) 
227        perturb();
228      int i;
229      if (scalingFlag_>0) {
230        for (i=0;i<numberRows_;i++) {
231          dual_[i] = saveDuals[i]/rowScale_[i];
232        }
233      } else {
234        memcpy(dual_,saveDuals,numberRows_*sizeof(double));
235      }
236      // now create my duals
237      for (i=0;i<numberRows_;i++) {
238        // slack
239        double value = dual_[i];
240        value += rowObjectiveWork_[i];
241        saveDuals[i+numberColumns_]=value;
242      }
243      ClpDisjointCopyN(objectiveWork_,numberColumns_,saveDuals);
244      transposeTimes(-1.0,dual_,saveDuals);
245      memcpy(dj_,saveDuals,(numberColumns_+numberRows_)*sizeof(double));
246      // set up possible ones
247      for (i=0;i<numberRows_+numberColumns_;i++)
248        clearPivoted(i);
249      int iRow;
250      for (iRow=0;iRow<numberRows_;iRow++) {
251        int iPivot=pivotVariable_[iRow];
252        if (fabs(saveDuals[iPivot])>dualTolerance_) {
253          if (getStatus(iPivot)!=isFree) 
254            setPivoted(iPivot);
255        }
256      }
257    } 
258
259    double objectiveChange;
260    numberFake_ =0; // Number of variables at fake bounds
261    changeBounds(true,NULL,objectiveChange);
262   
263    int lastCleaned=0; // last time objective or bounds cleaned up
264
265    if (!ifValuesPass) {
266      // Check optimal
267      if (!numberDualInfeasibilities_&&!numberPrimalInfeasibilities_)
268        problemStatus_=0;
269    }
270    if (problemStatus_<0&&perturbation_<100) {
271      perturb();
272      // Can't get here if values pass
273      gutsOfSolution(NULL,NULL);
274      if (handler_->logLevel()>2) {
275        handler_->message(CLP_SIMPLEX_STATUS,messages_)
276          <<numberIterations_<<objectiveValue();
277        handler_->printing(sumPrimalInfeasibilities_>0.0)
278          <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
279        handler_->printing(sumDualInfeasibilities_>0.0)
280          <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
281        handler_->printing(numberDualInfeasibilitiesWithoutFree_
282                           <numberDualInfeasibilities_)
283                             <<numberDualInfeasibilitiesWithoutFree_;
284        handler_->message()<<CoinMessageEol;
285      }
286    }
287       
288    // This says whether to restore things etc
289    int factorType=0;
290    /*
291      Status of problem:
292      0 - optimal
293      1 - infeasible
294      2 - unbounded
295      -1 - iterating
296      -2 - factorization wanted
297      -3 - redo checking without factorization
298      -4 - looks infeasible
299    */
300    while (problemStatus_<0) {
301      int iRow, iColumn;
302      // clear
303      for (iRow=0;iRow<4;iRow++) {
304        rowArray_[iRow]->clear();
305      }   
306     
307      for (iColumn=0;iColumn<2;iColumn++) {
308        columnArray_[iColumn]->clear();
309      }   
310     
311      // give matrix (and model costs and bounds a chance to be
312      // refreshed (normally null)
313      matrix_->refresh(this);
314      // If getting nowhere - why not give it a kick
315      // does not seem to work too well - do some more work
316      if (perturbation_<101&&numberIterations_>2*(numberRows_+numberColumns_)
317          &&initialStatus!=10) {
318        perturb();
319        // Can't get here if values pass
320        gutsOfSolution(NULL,NULL);
321        if (handler_->logLevel()>2) {
322          handler_->message(CLP_SIMPLEX_STATUS,messages_)
323            <<numberIterations_<<objectiveValue();
324          handler_->printing(sumPrimalInfeasibilities_>0.0)
325            <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
326          handler_->printing(sumDualInfeasibilities_>0.0)
327            <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
328          handler_->printing(numberDualInfeasibilitiesWithoutFree_
329                             <numberDualInfeasibilities_)
330                               <<numberDualInfeasibilitiesWithoutFree_;
331          handler_->message()<<CoinMessageEol;
332        }
333      }
334      // may factorize, checks if problem finished
335      statusOfProblemInDual(lastCleaned,factorType,progress_,saveDuals);
336      // If values pass then do easy ones on first time
337      if (ifValuesPass&&
338          progress_->lastIterationNumber(0)<0) {
339        doEasyOnesInValuesPass(saveDuals);
340      }
341     
342      // Say good factorization
343      factorType=1;
344      if (data.sparseThreshold_) {
345        // use default at present
346        factorization_->sparseThreshold(0);
347        factorization_->goSparse();
348      }
349
350      // exit if victory declared
351      if (problemStatus_>=0)
352        break;
353     
354      // Do iterations
355      whileIterating(saveDuals);
356    }
357  }
358
359  assert (problemStatus_||!sumPrimalInfeasibilities_);
360
361  // clean up
362  finish();
363  delete [] saveDuals;
364
365  // Restore any saved stuff
366  restoreData(data);
367  return problemStatus_;
368}
369/* Reasons to come out:
370   -1 iterations etc
371   -2 inaccuracy
372   -3 slight inaccuracy (and done iterations)
373   +0 looks optimal (might be unbounded - but we will investigate)
374   +1 looks infeasible
375   +3 max iterations
376 */
377int
378ClpSimplexDual::whileIterating(double * & givenDuals)
379{
380#ifdef CLP_DEBUG
381  int debugIteration=-1;
382#endif
383  {
384    int i;
385    for (i=0;i<4;i++) {
386      rowArray_[i]->clear();
387    }   
388    for (i=0;i<2;i++) {
389      columnArray_[i]->clear();
390    }   
391  }     
392  // if can't trust much and long way from optimal then relax
393  if (largestPrimalError_>10.0)
394    factorization_->relaxAccuracyCheck(min(1.0e2,largestPrimalError_/10.0));
395  else
396    factorization_->relaxAccuracyCheck(1.0);
397  // status stays at -1 while iterating, >=0 finished, -2 to invert
398  // status -3 to go to top without an invert
399  int returnCode = -1;
400
401  // compute average infeasibility for backward test
402  double averagePrimalInfeasibility = sumPrimalInfeasibilities_/
403    ((double ) (numberPrimalInfeasibilities_+1));
404
405  // Get dubious weights
406  factorization_->getWeights(rowArray_[0]->getIndices());
407  CoinBigIndex * dubiousWeights = matrix_->dubiousWeights(this,rowArray_[0]->getIndices());
408  // If values pass then get list of candidates
409  int * candidateList = NULL;
410  int numberCandidates = 0;
411#ifdef CLP_DEBUG
412  bool wasInValuesPass= (givenDuals!=NULL);
413#endif
414  int candidate=-1;
415  if (givenDuals) {
416    candidateList = new int[numberRows_];
417    // move reduced costs across
418    memcpy(dj_,givenDuals,(numberRows_+numberColumns_)*sizeof(double));
419    int iRow;
420    for (iRow=0;iRow<numberRows_;iRow++) {
421      int iPivot=pivotVariable_[iRow];
422      if (flagged(iPivot))
423        continue;
424      if (fabs(dj_[iPivot])>dualTolerance_) {
425        if (pivoted(iPivot)) 
426          candidateList[numberCandidates++]=iRow;
427      } else {
428        clearPivoted(iPivot);
429      }
430    }
431    // and set first candidate
432    if (!numberCandidates) {
433      delete [] candidateList;
434      delete [] givenDuals;
435      givenDuals=NULL;
436      candidateList=NULL;
437      int iRow;
438      for (iRow=0;iRow<numberRows_;iRow++) {
439        int iPivot=pivotVariable_[iRow];
440        clearPivoted(iPivot);
441      }
442    }
443  }
444  while (problemStatus_==-1) {
445#ifdef CLP_DEBUG
446    if (givenDuals) {
447      double value5=0.0;
448      int i;
449      for (i=0;i<numberRows_+numberColumns_;i++) {
450        if (dj_[i]<-1.0e-6)
451          if (upper_[i]<1.0e20)
452            value5 += dj_[i]*upper_[i];
453          else
454            printf("bad dj %g on %d with large upper status %d\n",
455                   dj_[i],i,status_[i]&7);
456        else if (dj_[i] >1.0e-6)
457          if (lower_[i]>-1.0e20)
458            value5 += dj_[i]*lower_[i];
459          else
460            printf("bad dj %g on %d with large lower status %d\n",
461                   dj_[i],i,status_[i]&7);
462      }
463      printf("Values objective Value %g\n",value5);
464    }
465    if ((handler_->logLevel()&32)&&wasInValuesPass) {
466      double value5=0.0;
467      int i;
468      for (i=0;i<numberRows_+numberColumns_;i++) {
469        if (dj_[i]<-1.0e-6)
470          if (upper_[i]<1.0e20)
471            value5 += dj_[i]*upper_[i];
472        else if (dj_[i] >1.0e-6)
473          if (lower_[i]>-1.0e20)
474            value5 += dj_[i]*lower_[i];
475      }
476      printf("Values objective Value %g\n",value5);
477      {
478        int i;
479        for (i=0;i<numberRows_+numberColumns_;i++) {
480          int iSequence = i;
481          double oldValue;
482         
483          switch(getStatus(iSequence)) {
484           
485          case basic:
486          case ClpSimplex::isFixed:
487            break;
488          case isFree:
489          case superBasic:
490            abort();
491            break;
492          case atUpperBound:
493            oldValue = dj_[iSequence];
494            //assert (oldValue<=tolerance);
495            assert (fabs(solution_[iSequence]-upper_[iSequence])<1.0e-7);
496            break;
497          case atLowerBound:
498            oldValue = dj_[iSequence];
499            //assert (oldValue>=-tolerance);
500            assert (fabs(solution_[iSequence]-lower_[iSequence])<1.0e-7);
501            break;
502          }
503        }
504      }
505    }
506#endif
507#ifdef CLP_DEBUG
508    {
509      int i;
510      for (i=0;i<4;i++) {
511        rowArray_[i]->checkClear();
512      }   
513      for (i=0;i<2;i++) {
514        columnArray_[i]->checkClear();
515      }   
516    }     
517#endif
518#if CLP_DEBUG>2
519    // very expensive
520    if (numberIterations_>0&&numberIterations_<-801) {
521      handler_->setLogLevel(63);
522      double saveValue = objectiveValue_;
523      double * saveRow1 = new double[numberRows_];
524      double * saveRow2 = new double[numberRows_];
525      memcpy(saveRow1,rowReducedCost_,numberRows_*sizeof(double));
526      memcpy(saveRow2,rowActivityWork_,numberRows_*sizeof(double));
527      double * saveColumn1 = new double[numberColumns_];
528      double * saveColumn2 = new double[numberColumns_];
529      memcpy(saveColumn1,reducedCostWork_,numberColumns_*sizeof(double));
530      memcpy(saveColumn2,columnActivityWork_,numberColumns_*sizeof(double));
531      gutsOfSolution(NULL,NULL);
532      printf("xxx %d old obj %g, recomputed %g, sum dual inf %g\n",
533             numberIterations_,
534             saveValue,objectiveValue_,sumDualInfeasibilities_);
535      if (saveValue>objectiveValue_+1.0e-2)
536        printf("**bad**\n");
537      memcpy(rowReducedCost_,saveRow1,numberRows_*sizeof(double));
538      memcpy(rowActivityWork_,saveRow2,numberRows_*sizeof(double));
539      memcpy(reducedCostWork_,saveColumn1,numberColumns_*sizeof(double));
540      memcpy(columnActivityWork_,saveColumn2,numberColumns_*sizeof(double));
541      delete [] saveRow1;
542      delete [] saveRow2;
543      delete [] saveColumn1;
544      delete [] saveColumn2;
545      objectiveValue_=saveValue;
546    }
547#endif
548#if 0
549    if (!factorization_->pivots()){
550      int iPivot;
551      double * array = rowArray_[3]->denseVector();
552      int i;
553      for (iPivot=0;iPivot<numberRows_;iPivot++) {
554        int iSequence = pivotVariable_[iPivot];
555        unpack(rowArray_[3],iSequence);
556        factorization_->updateColumn(rowArray_[2],rowArray_[3]);
557        assert (fabs(array[iPivot]-1.0)<1.0e-4);
558        array[iPivot]=0.0;
559        for (i=0;i<numberRows_;i++)
560          assert (fabs(array[i])<1.0e-4);
561        rowArray_[3]->clear();
562      }
563    }
564#endif
565#ifdef CLP_DEBUG
566    {
567      int iSequence, number=numberRows_+numberColumns_;
568      for (iSequence=0;iSequence<number;iSequence++) {
569        double lowerValue=lower_[iSequence];
570        double upperValue=upper_[iSequence];
571        double value=solution_[iSequence];
572        if(getStatus(iSequence)!=basic) {
573          assert(lowerValue>-1.0e20);
574          assert(upperValue<1.0e20);
575        }
576        switch(getStatus(iSequence)) {
577         
578        case basic:
579          break;
580        case isFree:
581        case superBasic:
582          break;
583        case atUpperBound:
584          assert (fabs(value-upperValue)<=primalTolerance_) ;
585          break;
586        case atLowerBound:
587        case ClpSimplex::isFixed:
588          assert (fabs(value-lowerValue)<=primalTolerance_) ;
589          break;
590        }
591      }
592    }
593    if(numberIterations_==debugIteration) {
594      printf("dodgy iteration coming up\n");
595    }
596#endif
597    // choose row to go out
598    // dualRow will go to virtual row pivot choice algorithm
599    // make sure values pass off if it should be
600    if (numberCandidates) 
601      candidate = candidateList[--numberCandidates];
602    else
603      candidate=-1;
604    dualRow(candidate);
605    if (pivotRow_>=0) {
606      // we found a pivot row
607      handler_->message(CLP_SIMPLEX_PIVOTROW,messages_)
608        <<pivotRow_
609        <<CoinMessageEol;
610      // check accuracy of weights
611      dualRowPivot_->checkAccuracy();
612      // Get good size for pivot
613      double acceptablePivot=1.0e-7;
614      if (factorization_->pivots()>10)
615        acceptablePivot=1.0e-5; // if we have iterated be more strict
616      else if (factorization_->pivots()>5)
617        acceptablePivot=1.0e-6; // if we have iterated be slightly more strict
618      // get sign for finding row of tableau
619      if (candidate<0) {
620        // normal iteration
621        // create as packed
622        double direction=directionOut_;
623        rowArray_[0]->createPacked(1,&pivotRow_,&direction);
624        factorization_->updateColumnTranspose(rowArray_[1],rowArray_[0]);
625        // put row of tableau in rowArray[0] and columnArray[0]
626        matrix_->transposeTimes(this,-1.0,
627                              rowArray_[0],rowArray_[3],columnArray_[0]);
628        // do ratio test for normal iteration
629        dualColumn(rowArray_[0],columnArray_[0],columnArray_[1],
630                 rowArray_[3],acceptablePivot,dubiousWeights);
631      } else {
632        // Make sure direction plausible
633        assert (upperOut_<1.0e50||lowerOut_>-1.0e50);
634        if (directionOut_<0&&fabs(valueOut_-upperOut_)>dualBound_+primalTolerance_) {
635          if (fabs(valueOut_-upperOut_)>fabs(valueOut_-lowerOut_))
636            directionOut_=1;
637        } else if (directionOut_>0&&fabs(valueOut_-upperOut_)<dualBound_+primalTolerance_) {
638          if (fabs(valueOut_-upperOut_)>fabs(valueOut_-lowerOut_))
639            directionOut_=-1;
640        }
641        double direction=directionOut_;
642        rowArray_[0]->createPacked(1,&pivotRow_,&direction);
643        factorization_->updateColumnTranspose(rowArray_[1],rowArray_[0]);
644        // put row of tableau in rowArray[0] and columnArray[0]
645        matrix_->transposeTimes(this,-1.0,
646                              rowArray_[0],rowArray_[3],columnArray_[0]);
647        acceptablePivot *= 10.0;
648        // do ratio test
649        checkPossibleValuesMove(rowArray_[0],columnArray_[0],
650                                            acceptablePivot,NULL);
651
652        // recompute true dualOut_
653        if (directionOut_<0) {
654          dualOut_ = valueOut_ - upperOut_;
655        } else {
656          dualOut_ = lowerOut_ - valueOut_;
657        }
658        // check what happened if was values pass
659        // may want to move part way i.e. movement
660        bool normalIteration = (sequenceIn_!=sequenceOut_);
661
662        clearPivoted(sequenceOut_);  // make sure won't be done again
663        // see if end of values pass
664        if (!numberCandidates) {
665          int iRow;
666          delete [] candidateList;
667          delete [] givenDuals;
668          candidate=-2; // -2 signals end
669          givenDuals=NULL;
670          handler_->message(CLP_END_VALUES_PASS,messages_)
671            <<numberIterations_;
672          candidateList=NULL;
673          for (iRow=0;iRow<numberRows_;iRow++) {
674            int iPivot=pivotVariable_[iRow];
675            //assert (fabs(dj_[iPivot]),1.0e-5);
676            clearPivoted(iPivot);
677          }
678        }
679        if (!normalIteration) {
680          //rowArray_[0]->cleanAndPackSafe(1.0e-60);
681          //columnArray_[0]->cleanAndPackSafe(1.0e-60);
682          updateDualsInValuesPass(rowArray_[0],columnArray_[0],theta_);
683          if (candidate==-2)
684            problemStatus_=-2;
685          continue; // skip rest of iteration
686        } else {
687          // recompute dualOut_
688          if (directionOut_<0) {
689            dualOut_ = valueOut_ - upperOut_;
690          } else {
691            dualOut_ = lowerOut_ - valueOut_;
692          }
693        }
694      }
695      if (sequenceIn_>=0) {
696        // normal iteration
697        // update the incoming column
698        double btranAlpha = -alpha_*directionOut_; // for check
699        unpackPacked(rowArray_[1]);
700        factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
701        // and update dual weights (can do in parallel - with extra array)
702        alpha_ = dualRowPivot_->updateWeights(rowArray_[0],
703                                              rowArray_[2],
704                                              rowArray_[1]);
705        // see if update stable
706#ifdef CLP_DEBUG
707        if ((handler_->logLevel()&32))
708          printf("btran alpha %g, ftran alpha %g\n",btranAlpha,alpha_);
709#endif
710        double checkValue=1.0e-7;
711        // if can't trust much and long way from optimal then relax
712        if (largestPrimalError_>10.0)
713          checkValue = min(1.0e-4,1.0e-8*largestPrimalError_);
714        if (fabs(btranAlpha)<1.0e-12||fabs(alpha_)<1.0e-12||
715            fabs(btranAlpha-alpha_)>checkValue*(1.0+fabs(alpha_))) {
716          handler_->message(CLP_DUAL_CHECK,messages_)
717            <<btranAlpha
718            <<alpha_
719            <<CoinMessageEol;
720          if (factorization_->pivots()) {
721            dualRowPivot_->unrollWeights();
722            problemStatus_=-2; // factorize now
723            rowArray_[0]->clear();
724            rowArray_[1]->clear();
725            columnArray_[0]->clear();
726            returnCode=-2;
727            break;
728          } else {
729            // take on more relaxed criterion
730            if (fabs(btranAlpha)<1.0e-12||fabs(alpha_)<1.0e-12||
731                fabs(btranAlpha-alpha_)>1.0e-4*(1.0+fabs(alpha_))) {
732              dualRowPivot_->unrollWeights();
733              // need to reject something
734              char x = isColumn(sequenceOut_) ? 'C' :'R';
735              handler_->message(CLP_SIMPLEX_FLAG,messages_)
736                <<x<<sequenceWithin(sequenceOut_)
737                <<CoinMessageEol;
738              setFlagged(sequenceOut_);
739              lastBadIteration_ = numberIterations_; // say be more cautious
740              rowArray_[0]->clear();
741              rowArray_[1]->clear();
742              columnArray_[0]->clear();
743              continue;
744            }
745          }
746        }
747        // update duals BEFORE replaceColumn so can do updateColumn
748        double objectiveChange=0.0;
749        // do duals first as variables may flip bounds
750        // rowArray_[0] and columnArray_[0] may have flips
751        // so use rowArray_[3] for work array from here on
752        int nswapped = 0;
753        //rowArray_[0]->cleanAndPackSafe(1.0e-60);
754        //columnArray_[0]->cleanAndPackSafe(1.0e-60);
755        if (candidate==-1) 
756          nswapped = updateDualsInDual(rowArray_[0],columnArray_[0],
757                                       rowArray_[2],theta_,
758                                       objectiveChange,false);
759        else 
760          updateDualsInValuesPass(rowArray_[0],columnArray_[0],theta_);
761
762        // which will change basic solution
763        if (nswapped) {
764          factorization_->updateColumn(rowArray_[3],rowArray_[2]);
765          dualRowPivot_->updatePrimalSolution(rowArray_[2],
766                                              1.0,objectiveChange);
767         
768          // recompute dualOut_
769          valueOut_ = solution_[sequenceOut_];
770          if (directionOut_<0) {
771            dualOut_ = valueOut_ - upperOut_;
772          } else {
773            dualOut_ = lowerOut_ - valueOut_;
774          }
775          if(dualOut_<-max(1.0e-12*averagePrimalInfeasibility,1.0e-8)
776             &&factorization_->pivots()&&
777             getStatus(sequenceIn_)!=isFree) {
778            // going backwards - factorize
779            dualRowPivot_->unrollWeights();
780            problemStatus_=-2; // factorize now
781            returnCode=-2;
782            break;
783          }
784        }
785        // amount primal will move
786        double movement = -dualOut_*directionOut_/alpha_;
787        // so objective should increase by fabs(dj)*movement
788        // but we already have objective change - so check will be good
789        if (objectiveChange+fabs(movement*dualIn_)<-1.0e-5) {
790#ifdef CLP_DEBUG
791          if (handler_->logLevel()&32)
792            printf("movement %g, swap change %g, rest %g  * %g\n",
793                   objectiveChange+fabs(movement*dualIn_),
794                   objectiveChange,movement,dualIn_);
795#endif
796          if(factorization_->pivots()>5) {
797            // going backwards - factorize
798            dualRowPivot_->unrollWeights();
799            problemStatus_=-2; // factorize now
800            returnCode=-2;
801            break;
802          }
803        }
804        assert(fabs(dualOut_)<1.0e50);
805        // if stable replace in basis
806        int updateStatus = factorization_->replaceColumn(rowArray_[2],
807                                                         pivotRow_,
808                                                         alpha_);
809        // if no pivots, bad update but reasonable alpha - take and invert
810        if (updateStatus==2&&
811                   !factorization_->pivots()&&fabs(alpha_)>1.0e-5)
812          updateStatus=4;
813        if (updateStatus==1||updateStatus==4) {
814          // slight error
815          if (factorization_->pivots()>5||updateStatus==4) {
816            problemStatus_=-2; // factorize now
817            returnCode=-3;
818          }
819        } else if (updateStatus==2) {
820          // major error
821          dualRowPivot_->unrollWeights();
822          // later we may need to unwind more e.g. fake bounds
823          if (factorization_->pivots()) {
824            problemStatus_=-2; // factorize now
825            returnCode=-2;
826            break;
827          } else {
828            // need to reject something
829            char x = isColumn(sequenceOut_) ? 'C' :'R';
830            handler_->message(CLP_SIMPLEX_FLAG,messages_)
831              <<x<<sequenceWithin(sequenceOut_)
832              <<CoinMessageEol;
833            setFlagged(sequenceOut_);
834            lastBadIteration_ = numberIterations_; // say be more cautious
835            rowArray_[0]->clear();
836            rowArray_[1]->clear();
837            columnArray_[0]->clear();
838            // make sure dual feasible
839            // look at all rows and columns
840            double objectiveChange=0.0;
841            updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[1],
842                              0.0,objectiveChange,true);
843            continue;
844          }
845        } else if (updateStatus==3) {
846          // out of memory
847          // increase space if not many iterations
848          if (factorization_->pivots()<
849              0.5*factorization_->maximumPivots()&&
850              factorization_->pivots()<200)
851            factorization_->areaFactor(
852                                       factorization_->areaFactor() * 1.1);
853          problemStatus_=-2; // factorize now
854        } else if (updateStatus==5) {
855          problemStatus_=-2; // factorize now
856        }
857        // update primal solution
858        if (theta_<0.0&&candidate==-1) {
859#ifdef CLP_DEBUG
860          if (handler_->logLevel()&32)
861            printf("negative theta %g\n",theta_);
862#endif
863          theta_=0.0;
864        }
865        // do actual flips
866        flipBounds(rowArray_[0],columnArray_[0],theta_);
867        //rowArray_[1]->expand();
868        dualRowPivot_->updatePrimalSolution(rowArray_[1],
869                                            movement,
870                                            objectiveChange);
871#ifdef CLP_DEBUG
872        double oldobj=objectiveValue_;
873#endif
874        // modify dualout
875        dualOut_ /= alpha_;
876        dualOut_ *= -directionOut_;
877        //setStatus(sequenceIn_,basic);
878        dj_[sequenceIn_]=0.0;
879        double oldValue=valueIn_;
880        if (directionIn_==-1) {
881          // as if from upper bound
882          valueIn_ = upperIn_+dualOut_;
883        } else {
884          // as if from lower bound
885          valueIn_ = lowerIn_+dualOut_;
886        }
887        objectiveChange += cost_[sequenceIn_]*(valueIn_-oldValue);
888        // outgoing
889        // set dj to zero unless values pass
890        if (directionOut_>0) {
891          valueOut_ = lowerOut_;
892          if (candidate==-1)
893            dj_[sequenceOut_] = theta_;
894        } else {
895          valueOut_ = upperOut_;
896          if (candidate==-1)
897            dj_[sequenceOut_] = -theta_;
898        }
899        solution_[sequenceOut_]=valueOut_;
900        int whatNext=housekeeping(objectiveChange);
901        // and set bounds correctly
902        originalBound(sequenceIn_); 
903        changeBound(sequenceOut_);
904#ifdef CLP_DEBUG
905        if (objectiveValue_<oldobj-1.0e-5&&(handler_->logLevel()&16))
906          printf("obj backwards %g %g\n",objectiveValue_,oldobj);
907#endif
908        if (whatNext==1||candidate==-2) {
909          problemStatus_ =-2; // refactorize
910        } else if (whatNext==2) {
911          // maximum iterations or equivalent
912          problemStatus_= 3;
913          returnCode=3;
914          break;
915        }
916      } else {
917        // no incoming column is valid
918#ifdef CLP_DEBUG
919        if (handler_->logLevel()&32)
920          printf("** no column pivot\n");
921#endif
922        if (factorization_->pivots()<5) {
923          problemStatus_=-4; //say looks infeasible
924          // create ray anyway
925          delete [] ray_;
926          ray_ = new double [ numberRows_];
927          rowArray_[0]->expand(); // in case packed
928          ClpDisjointCopyN(rowArray_[0]->denseVector(),numberRows_,ray_);
929        }
930        rowArray_[0]->clear();
931        columnArray_[0]->clear();
932        returnCode=1;
933        break;
934      }
935    } else {
936      // no pivot row
937#ifdef CLP_DEBUG
938      if (handler_->logLevel()&32)
939        printf("** no row pivot\n");
940#endif
941      if (!factorization_->pivots()) {
942        // may have crept through - so may be optimal
943        // check any flagged variables
944        int iRow;
945        for (iRow=0;iRow<numberRows_;iRow++) {
946          int iPivot=pivotVariable_[iRow];
947          if (flagged(iPivot))
948            break;
949        }
950        if (numberFake_||numberDualInfeasibilities_) {
951          // may be dual infeasible
952          problemStatus_=-5;
953        } else {
954          if (iRow<numberRows_) {
955#ifdef CLP_DEBUG
956            std::cerr<<"Flagged variables at end - infeasible?"<<std::endl;
957            printf("Probably infeasible - pivot was %g\n",alpha_);
958#endif
959            if (fabs(alpha_)<1.0e-4) {
960              problemStatus_=1;
961            } else {
962#ifdef CLP_DEBUG
963              abort();
964#endif
965            }
966          } else {
967            problemStatus_=0;
968          }
969        }
970      } else {
971        problemStatus_=-3;
972      }
973      returnCode=0;
974      break;
975    }
976  }
977  if (givenDuals) {
978    memcpy(givenDuals,dj_,(numberRows_+numberColumns_)*sizeof(double));
979    // get rid of any values pass array
980    delete [] candidateList;
981  }
982  delete [] dubiousWeights;
983  return returnCode;
984}
985/* The duals are updated by the given arrays.
986   Returns number of infeasibilities.
987   rowArray and columnarray will have flipped
988   The output vector has movement (row length array) */
989int
990ClpSimplexDual::updateDualsInDual(CoinIndexedVector * rowArray,
991                                  CoinIndexedVector * columnArray,
992                                  CoinIndexedVector * outputArray,
993                                  double theta,
994                                  double & objectiveChange,
995                                  bool fullRecompute)
996{
997 
998  outputArray->clear();
999 
1000 
1001  int numberInfeasibilities=0;
1002  int numberRowInfeasibilities=0;
1003 
1004  // use a tighter tolerance except for all being okay
1005  double tolerance = dualTolerance_;
1006 
1007  double changeObj=0.0;
1008
1009  // Coding is very similar but we can save a bit by splitting
1010  // Do rows
1011  if (!fullRecompute) {
1012    int i;
1013    double * reducedCost = djRegion(0);
1014    const double * lower = lowerRegion(0);
1015    const double * upper = upperRegion(0);
1016    const double * cost = costRegion(0);
1017    double * work;
1018    int number;
1019    int * which;
1020    assert(rowArray->packedMode());
1021    work = rowArray->denseVector();
1022    number = rowArray->getNumElements();
1023    which = rowArray->getIndices();
1024    for (i=0;i<number;i++) {
1025      int iSequence = which[i];
1026      double alphaI = work[i];
1027      work[i]=0.0;
1028     
1029      Status status = getStatus(iSequence+numberColumns_);
1030      // more likely to be at upper bound ?
1031      if (status==atUpperBound) {
1032
1033        double value = reducedCost[iSequence]-theta*alphaI;
1034        reducedCost[iSequence]=value;
1035
1036        if (value>tolerance) {
1037          // to lower bound (if swap)
1038          which[numberInfeasibilities++]=iSequence;
1039          double movement = lower[iSequence]-upper[iSequence];
1040#ifdef CLP_DEBUG
1041          if ((handler_->logLevel()&32))
1042            printf("%d %d, new dj %g, alpha %g, movement %g\n",
1043                   0,iSequence,value,alphaI,movement);
1044#endif
1045          changeObj += movement*cost[iSequence];
1046          outputArray->quickAdd(iSequence,-movement);
1047        }
1048      } else if (status==atLowerBound) {
1049
1050        double value = reducedCost[iSequence]-theta*alphaI;
1051        reducedCost[iSequence]=value;
1052
1053        if (value<-tolerance) {
1054          // to upper bound
1055          which[numberInfeasibilities++]=iSequence;
1056          double movement = upper[iSequence] - lower[iSequence];
1057#ifdef CLP_DEBUG
1058          if ((handler_->logLevel()&32))
1059            printf("%d %d, new dj %g, alpha %g, movement %g\n",
1060                   0,iSequence,value,alphaI,movement);
1061#endif
1062          changeObj += movement*cost[iSequence];
1063          outputArray->quickAdd(iSequence,-movement);
1064        }
1065      }
1066    }
1067  } else  {
1068    double * solution = solutionRegion(0);
1069    double * reducedCost = djRegion(0);
1070    const double * lower = lowerRegion(0);
1071    const double * upper = upperRegion(0);
1072    const double * cost = costRegion(0);
1073    int * which;
1074    which = rowArray->getIndices();
1075    int iSequence;
1076    for (iSequence=0;iSequence<numberRows_;iSequence++) {
1077      double value = reducedCost[iSequence];
1078     
1079      Status status = getStatus(iSequence+numberColumns_);
1080      // more likely to be at upper bound ?
1081      if (status==atUpperBound) {
1082        double movement=0.0;
1083
1084        if (value>tolerance) {
1085          // to lower bound (if swap)
1086          // put back alpha
1087          which[numberInfeasibilities++]=iSequence;
1088          movement = lower[iSequence]-upper[iSequence];
1089          changeObj += movement*cost[iSequence];
1090          outputArray->quickAdd(iSequence,-movement);
1091        } else if (value>-tolerance) {
1092          // at correct bound but may swap
1093          FakeBound bound = getFakeBound(iSequence+numberColumns_);
1094          if (bound==ClpSimplexDual::upperFake) {
1095            movement = lower[iSequence]-upper[iSequence];
1096            setStatus(iSequence+numberColumns_,atLowerBound);
1097            solution[iSequence] = lower[iSequence];
1098            changeObj += movement*cost[iSequence];
1099            numberFake_--;
1100            setFakeBound(iSequence+numberColumns_,noFake);
1101          }
1102        }
1103      } else if (status==atLowerBound) {
1104        double movement=0.0;
1105
1106        if (value<-tolerance) {
1107          // to upper bound
1108          // put back alpha
1109          which[numberInfeasibilities++]=iSequence;
1110          movement = upper[iSequence] - lower[iSequence];
1111          changeObj += movement*cost[iSequence];
1112          outputArray->quickAdd(iSequence,-movement);
1113        } else if (value<tolerance) {
1114          // at correct bound but may swap
1115          FakeBound bound = getFakeBound(iSequence+numberColumns_);
1116          if (bound==ClpSimplexDual::lowerFake) {
1117            movement = upper[iSequence]-lower[iSequence];
1118            setStatus(iSequence+numberColumns_,atUpperBound);
1119            solution[iSequence] = upper[iSequence];
1120            changeObj += movement*cost[iSequence];
1121            numberFake_--;
1122            setFakeBound(iSequence+numberColumns_,noFake);
1123          }
1124        }
1125      }
1126    }
1127  }
1128
1129  // Do columns
1130  if (!fullRecompute) {
1131    int i;
1132    double * reducedCost = djRegion(1);
1133    const double * lower = lowerRegion(1);
1134    const double * upper = upperRegion(1);
1135    const double * cost = costRegion(1);
1136    double * work;
1137    int number;
1138    int * which;
1139    // set number of infeasibilities in row array
1140    numberRowInfeasibilities=numberInfeasibilities;
1141    rowArray->setNumElements(numberInfeasibilities);
1142    numberInfeasibilities=0;
1143    work = columnArray->denseVector();
1144    number = columnArray->getNumElements();
1145    which = columnArray->getIndices();
1146   
1147    for (i=0;i<number;i++) {
1148      int iSequence = which[i];
1149      double alphaI = work[i];
1150      work[i]=0.0;
1151     
1152      Status status = getStatus(iSequence);
1153      if (status==atLowerBound) {
1154        double value = reducedCost[iSequence]-theta*alphaI;
1155        reducedCost[iSequence]=value;
1156        double movement=0.0;
1157
1158        if (value<-tolerance) {
1159          // to upper bound
1160          which[numberInfeasibilities++]=iSequence;
1161          movement = upper[iSequence] - lower[iSequence];
1162#ifdef CLP_DEBUG
1163          if ((handler_->logLevel()&32))
1164            printf("%d %d, new dj %g, alpha %g, movement %g\n",
1165                   1,iSequence,value,alphaI,movement);
1166#endif
1167          changeObj += movement*cost[iSequence];
1168          matrix_->add(this,outputArray,iSequence,movement);
1169        }
1170      } else if (status==atUpperBound) {
1171        double value = reducedCost[iSequence]-theta*alphaI;
1172        reducedCost[iSequence]=value;
1173        double movement=0.0;
1174
1175        if (value>tolerance) {
1176          // to lower bound (if swap)
1177          which[numberInfeasibilities++]=iSequence;
1178          movement = lower[iSequence]-upper[iSequence];
1179#ifdef CLP_DEBUG
1180          if ((handler_->logLevel()&32))
1181            printf("%d %d, new dj %g, alpha %g, movement %g\n",
1182                   1,iSequence,value,alphaI,movement);
1183#endif
1184          changeObj += movement*cost[iSequence];
1185          matrix_->add(this,outputArray,iSequence,movement);
1186        }
1187      } else if (status==isFree) {
1188        double value = reducedCost[iSequence]-theta*alphaI;
1189        reducedCost[iSequence]=value;
1190      }
1191    }
1192  } else {
1193    double * solution = solutionRegion(1);
1194    double * reducedCost = djRegion(1);
1195    const double * lower = lowerRegion(1);
1196    const double * upper = upperRegion(1);
1197    const double * cost = costRegion(1);
1198    int * which;
1199    // set number of infeasibilities in row array
1200    numberRowInfeasibilities=numberInfeasibilities;
1201    rowArray->setNumElements(numberInfeasibilities);
1202    numberInfeasibilities=0;
1203    which = columnArray->getIndices();
1204    int iSequence;
1205    for (iSequence=0;iSequence<numberColumns_;iSequence++) {
1206      double value = reducedCost[iSequence];
1207     
1208      Status status = getStatus(iSequence);
1209      if (status==atLowerBound) {
1210        double movement=0.0;
1211
1212        if (value<-tolerance) {
1213          // to upper bound
1214          // put back alpha
1215          which[numberInfeasibilities++]=iSequence;
1216          movement = upper[iSequence] - lower[iSequence];
1217          changeObj += movement*cost[iSequence];
1218          matrix_->add(this,outputArray,iSequence,movement);
1219        } else if (value<tolerance) {
1220          // at correct bound but may swap
1221          FakeBound bound = getFakeBound(iSequence);
1222          if (bound==ClpSimplexDual::lowerFake) {
1223            movement = upper[iSequence]-lower[iSequence];
1224            setStatus(iSequence,atUpperBound);
1225            solution[iSequence] = upper[iSequence];
1226            changeObj += movement*cost[iSequence];
1227            numberFake_--;
1228            setFakeBound(iSequence,noFake);
1229          }
1230        }
1231      } else if (status==atUpperBound) {
1232        double movement=0.0;
1233
1234        if (value>tolerance) {
1235          // to lower bound (if swap)
1236          // put back alpha
1237          which[numberInfeasibilities++]=iSequence;
1238          movement = lower[iSequence]-upper[iSequence];
1239          changeObj += movement*cost[iSequence];
1240          matrix_->add(this,outputArray,iSequence,movement);
1241        } else if (value>-tolerance) {
1242          // at correct bound but may swap
1243          FakeBound bound = getFakeBound(iSequence);
1244          if (bound==ClpSimplexDual::upperFake) {
1245            movement = lower[iSequence]-upper[iSequence];
1246            setStatus(iSequence,atLowerBound);
1247            solution[iSequence] = lower[iSequence];
1248            changeObj += movement*cost[iSequence];
1249            numberFake_--;
1250            setFakeBound(iSequence,noFake);
1251          }
1252        }
1253      }
1254    }
1255  }
1256#ifdef CLP_DEBUG
1257  if (fullRecompute&&numberFake_&&(handler_->logLevel()&16)!=0) 
1258    printf("%d fake after full update\n",numberFake_);
1259#endif
1260  // set number of infeasibilities
1261  columnArray->setNumElements(numberInfeasibilities);
1262  numberInfeasibilities += numberRowInfeasibilities;
1263  if (fullRecompute) {
1264    // do actual flips
1265    flipBounds(rowArray,columnArray,0.0);
1266  }
1267  objectiveChange += changeObj;
1268  return numberInfeasibilities;
1269}
1270void
1271ClpSimplexDual::updateDualsInValuesPass(CoinIndexedVector * rowArray,
1272                                  CoinIndexedVector * columnArray,
1273                                        double theta)
1274{
1275 
1276  // use a tighter tolerance except for all being okay
1277  double tolerance = dualTolerance_;
1278 
1279  // Coding is very similar but we can save a bit by splitting
1280  // Do rows
1281  {
1282    int i;
1283    double * reducedCost = djRegion(0);
1284    double * work;
1285    int number;
1286    int * which;
1287    work = rowArray->denseVector();
1288    number = rowArray->getNumElements();
1289    which = rowArray->getIndices();
1290    for (i=0;i<number;i++) {
1291      int iSequence = which[i];
1292      double alphaI = work[i];
1293      double value = reducedCost[iSequence]-theta*alphaI;
1294      work[i]=0.0;
1295      reducedCost[iSequence]=value;
1296     
1297      Status status = getStatus(iSequence+numberColumns_);
1298      // more likely to be at upper bound ?
1299      if (status==atUpperBound) {
1300
1301        if (value>tolerance) 
1302          reducedCost[iSequence]=0.0;
1303      } else if (status==atLowerBound) {
1304
1305        if (value<-tolerance) {
1306          reducedCost[iSequence]=0.0;
1307        }
1308      }
1309    }
1310  }
1311  rowArray->setNumElements(0);
1312
1313  // Do columns
1314  {
1315    int i;
1316    double * reducedCost = djRegion(1);
1317    double * work;
1318    int number;
1319    int * which;
1320    work = columnArray->denseVector();
1321    number = columnArray->getNumElements();
1322    which = columnArray->getIndices();
1323   
1324    for (i=0;i<number;i++) {
1325      int iSequence = which[i];
1326      double alphaI = work[i];
1327      double value = reducedCost[iSequence]-theta*alphaI;
1328      work[i]=0.0;
1329      reducedCost[iSequence]=value;
1330     
1331      Status status = getStatus(iSequence);
1332      if (status==atLowerBound) {
1333        if (value<-tolerance) 
1334          reducedCost[iSequence]=0.0;
1335      } else if (status==atUpperBound) {
1336        if (value>tolerance) 
1337          reducedCost[iSequence]=0.0;
1338      }
1339    }
1340  }
1341  columnArray->setNumElements(0);
1342}
1343/*
1344   Chooses dual pivot row
1345   Would be faster with separate region to scan
1346   and will have this (with square of infeasibility) when steepest
1347   For easy problems we can just choose one of the first rows we look at
1348*/
1349void
1350ClpSimplexDual::dualRow(int alreadyChosen)
1351{
1352  // get pivot row using whichever method it is
1353  int chosenRow=-1;
1354  if (alreadyChosen<0) {
1355    // first see if any free variables and put them in basis
1356    int nextFree = nextSuperBasic();
1357    //nextFree=-1; //off
1358    if (nextFree>=0) {
1359      // unpack vector and find a good pivot
1360      unpack(rowArray_[1],nextFree);
1361      factorization_->updateColumn(rowArray_[2],rowArray_[1]);
1362     
1363      double * work=rowArray_[1]->denseVector();
1364      int number=rowArray_[1]->getNumElements();
1365      int * which=rowArray_[1]->getIndices();
1366      double bestFeasibleAlpha=0.0;
1367      int bestFeasibleRow=-1;
1368      double bestInfeasibleAlpha=0.0;
1369      int bestInfeasibleRow=-1;
1370      int i;
1371     
1372      for (i=0;i<number;i++) {
1373        int iRow = which[i];
1374        double alpha = fabs(work[iRow]);
1375        if (alpha>1.0e-3) {
1376          int iSequence=pivotVariable_[iRow];
1377          double value = solution_[iSequence];
1378          double lower = lower_[iSequence];
1379          double upper = upper_[iSequence];
1380          double infeasibility=0.0;
1381          if (value>upper)
1382            infeasibility = value-upper;
1383          else if (value<lower)
1384            infeasibility = lower-value;
1385          if (infeasibility*alpha>bestInfeasibleAlpha&&alpha>1.0e-1) {
1386            if (!flagged(iSequence)) {
1387              bestInfeasibleAlpha = infeasibility*alpha;
1388              bestInfeasibleRow=iRow;
1389            }
1390          }
1391          if (alpha>bestFeasibleAlpha&&(lower>-1.0e20||upper<1.0e20)) {
1392            bestFeasibleAlpha = alpha;
1393            bestFeasibleRow=iRow;
1394          }
1395        }
1396      }
1397      if (bestInfeasibleRow>=0) 
1398        chosenRow = bestInfeasibleRow;
1399      else if (bestFeasibleAlpha>1.0e-2)
1400        chosenRow = bestFeasibleRow;
1401      if (chosenRow>=0)
1402        pivotRow_=chosenRow;
1403      rowArray_[1]->clear();
1404    } 
1405  } else {
1406    // in values pass
1407    chosenRow=alreadyChosen;
1408    pivotRow_=chosenRow;
1409  }
1410  if (chosenRow<0) 
1411    pivotRow_=dualRowPivot_->pivotRow();
1412
1413  if (pivotRow_>=0) {
1414    sequenceOut_ = pivotVariable_[pivotRow_];
1415    valueOut_ = solution_[sequenceOut_];
1416    lowerOut_ = lower_[sequenceOut_];
1417    upperOut_ = upper_[sequenceOut_];
1418    if (alreadyChosen<0) {
1419      // if we have problems we could try other way and hope we get a
1420      // zero pivot?
1421      if (valueOut_>upperOut_) {
1422        directionOut_ = -1;
1423        dualOut_ = valueOut_ - upperOut_;
1424      } else if (valueOut_<lowerOut_) {
1425        directionOut_ = 1;
1426        dualOut_ = lowerOut_ - valueOut_;
1427      } else {
1428        // odd (could be free) - it's feasible - go to nearest
1429        if (valueOut_-lowerOut_<upperOut_-valueOut_) {
1430          directionOut_ = 1;
1431          dualOut_ = lowerOut_ - valueOut_;
1432        } else {
1433          directionOut_ = -1;
1434          dualOut_ = valueOut_ - upperOut_;
1435        }
1436      }
1437#ifdef CLP_DEBUG
1438      assert(dualOut_>=0.0);
1439#endif
1440    } else {
1441      // in values pass so just use sign of dj
1442      // We don't want to go through any barriers so set dualOut low
1443      // free variables will never be here
1444      dualOut_ = 1.0e-6;
1445      if (dj_[sequenceOut_]>0.0) {
1446        // this will give a -1 in pivot row (as slacks are -1.0)
1447        directionOut_ = 1;
1448      } else {
1449        directionOut_ = -1;
1450      }
1451    }
1452  }
1453  return ;
1454}
1455// Checks if any fake bounds active - if so returns number and modifies
1456// dualBound_ and everything.
1457// Free variables will be left as free
1458// Returns number of bounds changed if >=0
1459// Returns -1 if not initialize and no effect
1460// Fills in changeVector which can be used to see if unbounded
1461// and cost of change vector
1462int
1463ClpSimplexDual::changeBounds(bool initialize,
1464                                 CoinIndexedVector * outputArray,
1465                                 double & changeCost)
1466{ 
1467  numberFake_ = 0;
1468  if (!initialize) {
1469    int numberInfeasibilities;
1470    double newBound;
1471    newBound = 5.0*dualBound_;
1472    numberInfeasibilities=0;
1473    changeCost=0.0;
1474    // put back original bounds and then check
1475    createRim(3);
1476    int iSequence;
1477    // bounds will get bigger - just look at ones at bounds
1478    for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
1479      double lowerValue=lower_[iSequence];
1480      double upperValue=upper_[iSequence];
1481      double value=solution_[iSequence];
1482      setFakeBound(iSequence,ClpSimplexDual::noFake);
1483      switch(getStatus(iSequence)) {
1484
1485      case basic:
1486      case ClpSimplex::isFixed:
1487        break;
1488      case isFree:
1489      case superBasic:
1490        break;
1491      case atUpperBound:
1492        if (fabs(value-upperValue)>primalTolerance_) 
1493          numberInfeasibilities++;
1494        break;
1495      case atLowerBound:
1496        if (fabs(value-lowerValue)>primalTolerance_) 
1497          numberInfeasibilities++;
1498        break;
1499      }
1500    }
1501    // If dual infeasible then carry on
1502    if (numberInfeasibilities) {
1503      handler_->message(CLP_DUAL_CHECKB,messages_)
1504        <<newBound
1505        <<CoinMessageEol;
1506      int iSequence;
1507      for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
1508        double lowerValue=lower_[iSequence];
1509        double upperValue=upper_[iSequence];
1510        double newLowerValue;
1511        double newUpperValue;
1512        Status status = getStatus(iSequence);
1513        if (status==atUpperBound||
1514            status==atLowerBound) {
1515          double value = solution_[iSequence];
1516          if (value-lowerValue<=upperValue-value) {
1517            newLowerValue = max(lowerValue,value-0.666667*newBound);
1518            newUpperValue = min(upperValue,newLowerValue+newBound);
1519          } else {
1520            newUpperValue = min(upperValue,value+0.666667*newBound);
1521            newLowerValue = max(lowerValue,newUpperValue-newBound);
1522          }
1523          lower_[iSequence]=newLowerValue;
1524          upper_[iSequence]=newUpperValue;
1525          if (newLowerValue > lowerValue) {
1526            if (newUpperValue < upperValue) {
1527              setFakeBound(iSequence,ClpSimplexDual::bothFake);
1528              numberFake_++;
1529            } else {
1530              setFakeBound(iSequence,ClpSimplexDual::lowerFake);
1531              numberFake_++;
1532            }
1533          } else {
1534            if (newUpperValue < upperValue) {
1535              setFakeBound(iSequence,ClpSimplexDual::upperFake);
1536              numberFake_++;
1537            }
1538          }
1539          if (status==atUpperBound)
1540            solution_[iSequence] = newUpperValue;
1541          else
1542            solution_[iSequence] = newLowerValue;
1543          double movement = solution_[iSequence] - value;
1544          if (movement&&outputArray) {
1545            if (iSequence>=numberColumns_) {
1546              outputArray->quickAdd(iSequence,-movement);
1547              changeCost += movement*cost_[iSequence];
1548            } else {
1549              matrix_->add(this,outputArray,iSequence,movement);
1550              changeCost += movement*cost_[iSequence];
1551            }
1552          }
1553        }
1554      }
1555      dualBound_ = newBound;
1556    } else {
1557      numberInfeasibilities=-1;
1558    }
1559    return numberInfeasibilities;
1560  } else {
1561    int iSequence;
1562     
1563    for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
1564      Status status = getStatus(iSequence);
1565      if (status==atUpperBound||
1566          status==atLowerBound) {
1567        double lowerValue=lower_[iSequence];
1568        double upperValue=upper_[iSequence];
1569        double value = solution_[iSequence];
1570        if (lowerValue>-largeValue_||upperValue<largeValue_) {
1571          if (lowerValue-value>-0.5*dualBound_||
1572              upperValue-value<0.5*dualBound_) {
1573            if (fabs(lowerValue-value)<=fabs(upperValue-value)) {
1574              if (upperValue > lowerValue + dualBound_) {
1575                upper_[iSequence]=lowerValue+dualBound_;
1576                setFakeBound(iSequence,ClpSimplexDual::upperFake);
1577                numberFake_++;
1578              }
1579            } else {
1580              if (lowerValue < upperValue - dualBound_) {
1581                lower_[iSequence]=upperValue-dualBound_;
1582                setFakeBound(iSequence,ClpSimplexDual::lowerFake);
1583                numberFake_++;
1584              }
1585            }
1586          } else {
1587            lower_[iSequence]=-0.5*dualBound_;
1588            upper_[iSequence]= 0.5*dualBound_;
1589            setFakeBound(iSequence,ClpSimplexDual::bothFake);
1590            numberFake_++;
1591          }
1592        } else {
1593          // set non basic free variables to fake bounds
1594          // I don't think we should ever get here
1595          assert(!("should not be here"));
1596          lower_[iSequence]=-0.5*dualBound_;
1597          upper_[iSequence]= 0.5*dualBound_;
1598          setFakeBound(iSequence,ClpSimplexDual::bothFake);
1599          numberFake_++;
1600          setStatus(iSequence,atUpperBound);
1601          solution_[iSequence]=0.5*dualBound_;
1602        }
1603      }
1604    }
1605
1606    return 1;
1607  }
1608}
1609/*
1610   Row array has row part of pivot row (as duals so sign may be switched)
1611   Column array has column part.
1612   This chooses pivot column.
1613   Spare array will be needed when we start getting clever.
1614   We will check for basic so spare array will never overflow.
1615   If necessary will modify costs
1616*/
1617void
1618ClpSimplexDual::dualColumn(CoinIndexedVector * rowArray,
1619                           CoinIndexedVector * columnArray,
1620                           CoinIndexedVector * spareArray,
1621                           CoinIndexedVector * spareArray2,
1622                           double acceptablePivot,
1623                           CoinBigIndex * dubiousWeights)
1624{
1625  double * work;
1626  int number;
1627  int * which;
1628  double * reducedCost;
1629  int iSection;
1630
1631  sequenceIn_=-1;
1632  int numberPossiblySwapped=0;
1633  int numberRemaining=0;
1634 
1635  double totalThru=0.0; // for when variables flip
1636
1637  double bestEverPivot=acceptablePivot;
1638  int lastSequence = -1;
1639  double lastPivot=0.0;
1640  double upperTheta;
1641  double newTolerance = dualTolerance_;
1642  //newTolerance = dualTolerance_+1.0e-6*dblParam_[ClpDualTolerance];
1643  // will we need to increase tolerance
1644  bool thisIncrease=false;
1645  // If we think we need to modify costs (not if something from broad sweep)
1646  bool modifyCosts=false;
1647  // Increase in objective due to swapping bounds (may be negative)
1648  double increaseInObjective=0.0;
1649
1650  // use spareArrays to put ones looked at in
1651  // we are going to flip flop between
1652  int iFlip = 0;
1653  // Possible list of pivots
1654  int interesting[2];
1655  // where possible swapped ones are
1656  int swapped[2];
1657  // for zeroing out arrays after
1658  int marker[2][2];
1659  // pivot elements
1660  double * array[2], * spare, * spare2;
1661  // indices
1662  int * indices[2], * index, * index2;
1663  spareArray->clear();
1664  spareArray2->clear();
1665  array[0] = spareArray->denseVector();
1666  indices[0] = spareArray->getIndices();
1667  spare = array[0];
1668  index = indices[0];
1669  array[1] = spareArray2->denseVector();
1670  indices[1] = spareArray2->getIndices();
1671  int i;
1672  const double * lower;
1673  const double * upper;
1674
1675  // initialize lists
1676  for (i=0;i<2;i++) {
1677    interesting[i]=0;
1678    swapped[i]=numberColumns_;
1679    marker[i][0]=0;
1680    marker[i][1]=numberColumns_;
1681  }
1682
1683  /*
1684    First we get a list of possible pivots.  We can also see if the
1685    problem looks infeasible or whether we want to pivot in free variable.
1686    This may make objective go backwards but can only happen a finite
1687    number of times and I do want free variables basic.
1688
1689    Then we flip back and forth.  At the start of each iteration
1690    interesting[iFlip] should have possible candidates and swapped[iFlip]
1691    will have pivots if we decide to take a previous pivot.
1692    At end of each iteration interesting[1-iFlip] should have
1693    candidates if we go through this theta and swapped[1-iFlip]
1694    pivots if we don't go through.
1695
1696    At first we increase theta and see what happens.  We start
1697    theta at a reasonable guess.  If in right area then we do bit by bit.
1698
1699   */
1700
1701  // do first pass to get possibles
1702  // We can also see if infeasible or pivoting on free
1703  double tentativeTheta = 1.0e22;
1704  upperTheta = 1.0e31;
1705  double freePivot = acceptablePivot;
1706
1707  for (iSection=0;iSection<2;iSection++) {
1708
1709    int addSequence;
1710
1711    if (!iSection) {
1712      lower = rowLowerWork_;
1713      upper = rowUpperWork_;
1714      work = rowArray->denseVector();
1715      number = rowArray->getNumElements();
1716      which = rowArray->getIndices();
1717      reducedCost = rowReducedCost_;
1718      addSequence = numberColumns_;
1719    } else {
1720      lower = columnLowerWork_;
1721      upper = columnUpperWork_;
1722      work = columnArray->denseVector();
1723      number = columnArray->getNumElements();
1724      which = columnArray->getIndices();
1725      reducedCost = reducedCostWork_;
1726      addSequence = 0;
1727    }
1728   
1729    for (i=0;i<number;i++) {
1730      int iSequence = which[i];
1731      double alpha;
1732      double oldValue;
1733      double value;
1734      bool keep;
1735
1736      switch(getStatus(iSequence+addSequence)) {
1737         
1738      case basic:
1739      case ClpSimplex::isFixed:
1740        break;
1741      case isFree:
1742      case superBasic:
1743        alpha = work[i];
1744        oldValue = reducedCost[iSequence];
1745        if (oldValue>dualTolerance_) {
1746          keep = true;
1747        } else if (oldValue<-dualTolerance_) {
1748          keep = true;
1749        } else {
1750          if (fabs(alpha)>10.0*acceptablePivot) 
1751            keep = true;
1752          else
1753            keep=false;
1754        }
1755        if (keep) {
1756          // free - choose largest
1757          if (fabs(alpha)>freePivot) {
1758            freePivot=fabs(alpha);
1759            sequenceIn_ = iSequence + addSequence;
1760            theta_=oldValue/alpha;
1761            alpha_=alpha;
1762          }
1763        }
1764        break;
1765      case atUpperBound:
1766        alpha = work[i];
1767        oldValue = reducedCost[iSequence];
1768        value = oldValue-tentativeTheta*alpha;
1769        assert (oldValue<=dualTolerance_*1.0001);
1770        if (value>newTolerance) {
1771          value = oldValue-upperTheta*alpha;
1772          if (value>newTolerance && -alpha>=acceptablePivot) 
1773            upperTheta = (oldValue-newTolerance)/alpha;
1774          // add to list
1775          spare[numberRemaining]=alpha;
1776          index[numberRemaining++]=iSequence+addSequence;
1777        }
1778        break;
1779      case atLowerBound:
1780        alpha = work[i];
1781        oldValue = reducedCost[iSequence];
1782        value = oldValue-tentativeTheta*alpha;
1783        assert (oldValue>=-dualTolerance_*1.0001);
1784        if (value<-newTolerance) {
1785          value = oldValue-upperTheta*alpha;
1786          if (value<-newTolerance && alpha>=acceptablePivot) 
1787            upperTheta = (oldValue+newTolerance)/alpha;
1788          // add to list
1789          spare[numberRemaining]=alpha;
1790          index[numberRemaining++]=iSequence+addSequence;
1791        }
1792        break;
1793      }
1794    }
1795  }
1796  interesting[0]=numberRemaining;
1797  marker[0][0] = numberRemaining;
1798
1799  if (!numberRemaining&&sequenceIn_<0)
1800    return ; // Looks infeasible
1801
1802  if (sequenceIn_>=0) {
1803    // free variable - always choose
1804  } else {
1805
1806    theta_=1.0e50;
1807    // now flip flop between spare arrays until reasonable theta
1808    tentativeTheta = max(10.0*upperTheta,1.0e-7);
1809   
1810    // loops increasing tentative theta until can't go through
1811   
1812    while (tentativeTheta < 1.0e22) {
1813      double thruThis = 0.0;
1814     
1815      double bestPivot=acceptablePivot;
1816      int bestSequence=-1;
1817     
1818      numberPossiblySwapped = numberColumns_;
1819      numberRemaining = 0;
1820
1821      upperTheta = 1.0e50;
1822
1823      spare = array[iFlip];
1824      index = indices[iFlip];
1825      spare2 = array[1-iFlip];
1826      index2 = indices[1-iFlip];
1827
1828      // try 3 different ways
1829      // 1 bias increase by ones with slightly wrong djs
1830      // 2 bias by all
1831      // 3 bias by all - tolerance (doesn't seem very good)
1832#define TRYBIAS 1
1833
1834
1835      double increaseInThis=0.0; //objective increase in this loop
1836     
1837      for (i=0;i<interesting[iFlip];i++) {
1838        int iSequence = index[i];
1839        double alpha = spare[i];
1840        double oldValue = dj_[iSequence];
1841        double value = oldValue-tentativeTheta*alpha;
1842
1843        if (alpha < 0.0) {
1844          //at upper bound
1845          if (value>newTolerance) {
1846            double range = upper_[iSequence] - lower_[iSequence];
1847            thruThis -= range*alpha;
1848#if TRYBIAS==1
1849            if (oldValue>0.0)
1850              increaseInThis -= oldValue*range;
1851#elif TRYBIAS==2
1852            increaseInThis -= oldValue*range;
1853#else
1854            increaseInThis -= (oldValue+dualTolerance_)*range;
1855#endif
1856            // goes on swapped list (also means candidates if too many)
1857            spare2[--numberPossiblySwapped]=alpha;
1858            index2[numberPossiblySwapped]=iSequence;
1859            if (fabs(alpha)>bestPivot) {
1860              bestPivot=fabs(alpha);
1861              bestSequence=numberPossiblySwapped;
1862            }
1863          } else {
1864            value = oldValue-upperTheta*alpha;
1865            if (value>newTolerance && -alpha>=acceptablePivot) 
1866              upperTheta = (oldValue-newTolerance)/alpha;
1867            spare2[numberRemaining]=alpha;
1868            index2[numberRemaining++]=iSequence;
1869          }
1870        } else {
1871          // at lower bound
1872          if (value<-newTolerance) {
1873            double range = upper_[iSequence] - lower_[iSequence];
1874            thruThis += range*alpha;
1875            //?? is this correct - and should we look at good ones
1876#if TRYBIAS==1
1877            if (oldValue<0.0)
1878              increaseInThis += oldValue*range;
1879#elif TRYBIAS==2
1880            increaseInThis += oldValue*range;
1881#else
1882            increaseInThis += (oldValue-dualTolerance_)*range;
1883#endif
1884            // goes on swapped list (also means candidates if too many)
1885            spare2[--numberPossiblySwapped]=alpha;
1886            index2[numberPossiblySwapped]=iSequence;
1887            if (fabs(alpha)>bestPivot) {
1888              bestPivot=fabs(alpha);
1889              bestSequence=numberPossiblySwapped;
1890            }
1891          } else {
1892            value = oldValue-upperTheta*alpha;
1893            if (value<-newTolerance && alpha>=acceptablePivot) 
1894              upperTheta = (oldValue+newTolerance)/alpha;
1895            spare2[numberRemaining]=alpha;
1896            index2[numberRemaining++]=iSequence;
1897          }
1898        }
1899      }
1900      swapped[1-iFlip]=numberPossiblySwapped;
1901      interesting[1-iFlip]=numberRemaining;
1902      marker[1-iFlip][0]= max(marker[1-iFlip][0],numberRemaining);
1903      marker[1-iFlip][1]= min(marker[1-iFlip][1],numberPossiblySwapped);
1904     
1905      if (totalThru+thruThis>=fabs(dualOut_)||
1906          increaseInObjective+increaseInThis<0.0) {
1907        // We should be pivoting in this batch
1908        // so compress down to this lot
1909        numberRemaining=0;
1910        for (i=numberColumns_-1;i>=swapped[1-iFlip];i--) {
1911          spare[numberRemaining]=spare2[i];
1912          index[numberRemaining++]=index2[i];
1913        }
1914        interesting[iFlip]=numberRemaining;
1915        int iTry;
1916#define MAXTRY 100
1917        // first get ratio with tolerance
1918        for (iTry=0;iTry<MAXTRY;iTry++) {
1919         
1920          upperTheta=1.0e50;
1921          numberPossiblySwapped = numberColumns_;
1922          numberRemaining = 0;
1923
1924          increaseInThis=0.0; //objective increase in this loop
1925
1926          thruThis=0.0;
1927
1928          spare = array[iFlip];
1929          index = indices[iFlip];
1930          spare2 = array[1-iFlip];
1931          index2 = indices[1-iFlip];
1932     
1933          for (i=0;i<interesting[iFlip];i++) {
1934            int iSequence=index[i];
1935            double alpha=spare[i];
1936            double oldValue = dj_[iSequence];
1937            double value = oldValue-upperTheta*alpha;
1938           
1939            if (alpha < 0.0) {
1940              //at upper bound
1941              if (value>newTolerance) {
1942                if (-alpha>=acceptablePivot) {
1943                  upperTheta = (oldValue-newTolerance)/alpha;
1944                  // recompute value and make sure works
1945                  value = oldValue-upperTheta*alpha;
1946                  if (value<0)
1947                    upperTheta *= 1.0 +1.0e-11; // must be large
1948                }
1949              }
1950            } else {
1951              // at lower bound
1952              if (value<-newTolerance) {
1953                if (alpha>=acceptablePivot) {
1954                  upperTheta = (oldValue+newTolerance)/alpha;
1955                  // recompute value and make sure works
1956                  value = oldValue-upperTheta*alpha;
1957                  if (value>0)
1958                    upperTheta *= 1.0 +1.0e-11; // must be large
1959                }
1960              }
1961            }
1962          }
1963          bestPivot=acceptablePivot;
1964          sequenceIn_=-1;
1965          double bestWeight=COIN_DBL_MAX;
1966          double largestPivot=acceptablePivot;
1967          // now choose largest and sum all ones which will go through
1968          //printf("XX it %d number %d\n",numberIterations_,interesting[iFlip]);
1969          for (i=0;i<interesting[iFlip];i++) {
1970            int iSequence=index[i];
1971            double alpha=spare[i];
1972            double value = dj_[iSequence]-upperTheta*alpha;
1973            double badDj=0.0;
1974           
1975            bool addToSwapped=false;
1976           
1977            if (alpha < 0.0) {
1978              //at upper bound
1979              if (value>=0.0) { 
1980                addToSwapped=true;
1981#if TRYBIAS==1
1982                badDj = -max(dj_[iSequence],0.0);
1983#elif TRYBIAS==2
1984                badDj = -dj_[iSequence];
1985#else
1986                badDj = -dj_[iSequence]-dualTolerance_;
1987#endif
1988              }
1989            } else {
1990              // at lower bound
1991              if (value<=0.0) {
1992                addToSwapped=true;
1993#if TRYBIAS==1
1994                badDj = min(dj_[iSequence],0.0);
1995#elif TRYBIAS==2
1996                badDj = dj_[iSequence];
1997#else
1998                badDj = dj_[iSequence]-dualTolerance_;
1999#endif
2000              }
2001            }
2002            if (!addToSwapped) {
2003              // add to list of remaining
2004              spare2[numberRemaining]=alpha;
2005              index2[numberRemaining++]=iSequence;
2006            } else {
2007              // add to list of swapped
2008              spare2[--numberPossiblySwapped]=alpha;
2009              index2[numberPossiblySwapped]=iSequence;
2010              // select if largest pivot
2011              bool take=false;
2012              double absAlpha = fabs(alpha);
2013              // User could do anything to break ties here
2014              double weight;
2015              if (dubiousWeights)
2016                weight=dubiousWeights[iSequence];
2017              else
2018                weight=1.0;
2019              weight += CoinDrand48()*1.0e-2;
2020              if (absAlpha>2.0*bestPivot) {
2021                take=true;
2022              } else if (absAlpha>0.5*largestPivot) {
2023                // could multiply absAlpha and weight
2024                if (weight*bestPivot<bestWeight*absAlpha)
2025                  take=true;
2026              }
2027              if (take) {
2028                sequenceIn_ = numberPossiblySwapped;
2029                bestPivot =  absAlpha;
2030                theta_ = dj_[iSequence]/alpha;
2031                largestPivot = max(largestPivot,bestPivot);
2032                bestWeight = weight;
2033                //printf(" taken seq %d alpha %g weight %d\n",
2034                //   iSequence,absAlpha,dubiousWeights[iSequence]);
2035              } else {
2036                //printf(" not taken seq %d alpha %g weight %d\n",
2037                //   iSequence,absAlpha,dubiousWeights[iSequence]);
2038              }
2039              double range = upper[iSequence] - lower[iSequence];
2040              thruThis += range*fabs(alpha);
2041              increaseInThis += badDj*range;
2042            }
2043          }
2044          swapped[1-iFlip]=numberPossiblySwapped;
2045          interesting[1-iFlip]=numberRemaining;
2046          marker[1-iFlip][0]= max(marker[1-iFlip][0],numberRemaining);
2047          marker[1-iFlip][1]= min(marker[1-iFlip][1],numberPossiblySwapped);
2048          // If we stop now this will be increase in objective (I think)
2049          double increase = (fabs(dualOut_)-totalThru)*theta_;
2050          increase += increaseInObjective;
2051          if (theta_<0.0)
2052            thruThis += fabs(dualOut_); // force using this one
2053          if (increaseInObjective<0.0&&increase<0.0&&lastSequence>=0) {
2054            // back
2055            // We may need to be more careful - we could do by
2056            // switch so we always do fine grained?
2057            bestPivot=0.0;
2058          } else {
2059            // add in
2060            totalThru += thruThis;
2061            increaseInObjective += increaseInThis;
2062          }
2063          if (bestPivot<0.1*bestEverPivot&&
2064              bestEverPivot>1.0e-6&&
2065              (bestPivot<1.0e-3||totalThru*2.0>fabs(dualOut_))) {
2066            // back to previous one
2067            sequenceIn_=lastSequence;
2068            // swap regions
2069            iFlip = 1-iFlip;
2070            break;
2071          } else if (sequenceIn_==-1&&upperTheta>largeValue_) {
2072            if (lastPivot>acceptablePivot) {
2073              // back to previous one
2074              sequenceIn_=lastSequence;
2075              // swap regions
2076              iFlip = 1-iFlip;
2077            } else {
2078              // can only get here if all pivots too small
2079            }
2080            break;
2081          } else if (totalThru>=fabs(dualOut_)) {
2082            modifyCosts=true; // fine grain - we can modify costs
2083            break; // no point trying another loop
2084          } else {
2085            lastSequence=sequenceIn_;
2086            if (bestPivot>bestEverPivot)
2087              bestEverPivot=bestPivot;
2088            iFlip = 1 -iFlip;
2089            modifyCosts=true; // fine grain - we can modify costs
2090          }
2091        }
2092        if (iTry==MAXTRY)
2093          iFlip = 1-iFlip; // flip back
2094        break;
2095      } else {
2096        // skip this lot
2097        if (bestPivot>1.0e-3||bestPivot>bestEverPivot) {
2098          bestEverPivot=bestPivot;
2099          lastSequence=bestSequence;
2100        } else {
2101          // keep old swapped
2102          memcpy(array[1-iFlip]+swapped[iFlip],
2103                 array[iFlip]+swapped[iFlip],
2104                 (numberColumns_-swapped[iFlip])*sizeof(double));
2105          memcpy(indices[1-iFlip]+swapped[iFlip],
2106                 indices[iFlip]+swapped[iFlip],
2107                 (numberColumns_-swapped[iFlip])*sizeof(int));
2108          marker[1-iFlip][1] = min(marker[1-iFlip][1],swapped[iFlip]);
2109          swapped[1-iFlip]=swapped[iFlip];
2110        }
2111        increaseInObjective += increaseInThis;
2112        iFlip = 1 - iFlip; // swap regions
2113        tentativeTheta = 2.0*upperTheta;
2114        totalThru += thruThis;
2115      }
2116    }
2117   
2118    // can get here without sequenceIn_ set but with lastSequence
2119    if (sequenceIn_<0&&lastSequence>=0) {
2120      // back to previous one
2121      sequenceIn_=lastSequence;
2122      // swap regions
2123      iFlip = 1-iFlip;
2124    }
2125   
2126#define MINIMUMTHETA 1.0e-12
2127    // Movement should be minimum for anti-degeneracy - unless
2128    // fixed variable out
2129    double minimumTheta;
2130    if (upperOut_>lowerOut_)
2131      minimumTheta=MINIMUMTHETA;
2132    else
2133      minimumTheta=0.0;
2134    if (sequenceIn_>=0) {
2135      // at this stage sequenceIn_ is just pointer into index array
2136      // flip just so we can use iFlip
2137      iFlip = 1 -iFlip;
2138      spare = array[iFlip];
2139      index = indices[iFlip];
2140      double oldValue;
2141      alpha_ = spare[sequenceIn_];
2142      sequenceIn_ = indices[iFlip][sequenceIn_];
2143      oldValue = dj_[sequenceIn_];
2144      theta_ = oldValue/alpha_;
2145      if (theta_<minimumTheta) {
2146        // can't pivot to zero
2147        if (oldValue-minimumTheta*alpha_>=-dualTolerance_) {
2148          theta_=minimumTheta;
2149        } else if (oldValue-minimumTheta*alpha_>=-newTolerance) {
2150          theta_=minimumTheta;
2151          thisIncrease=true;
2152        } else {
2153          theta_=(oldValue+newTolerance)/alpha_;
2154          assert(theta_>=0.0);
2155          thisIncrease=true;
2156        } 
2157      }
2158      // may need to adjust costs so all dual feasible AND pivoted is exactly 0
2159      //int costOffset = numberRows_+numberColumns_;
2160      if (modifyCosts) {
2161        int i;
2162        for (i=numberColumns_-1;i>=swapped[iFlip];i--) {
2163          int iSequence=index[i];
2164          double alpha=spare[i];
2165          double value = dj_[iSequence]-theta_*alpha;
2166           
2167          // can't be free here
2168         
2169          if (alpha < 0.0) {
2170            //at upper bound
2171            if (value>dualTolerance_) {
2172              thisIncrease=true;
2173#define MODIFYCOST 2
2174#if MODIFYCOST
2175              // modify cost to hit new tolerance
2176              double modification = alpha*theta_-dj_[iSequence]
2177                +newTolerance;
2178              //modification = min(modification,dualTolerance_);
2179              //assert (fabs(modification)<1.0e-7);
2180              dj_[iSequence] += modification;
2181              cost_[iSequence] +=  modification;
2182              //cost_[iSequence+costOffset] += modification; // save change
2183#endif
2184            }
2185          } else {
2186            // at lower bound
2187            if (-value>dualTolerance_) {
2188              thisIncrease=true;
2189#if MODIFYCOST
2190              // modify cost to hit new tolerance
2191              double modification = alpha*theta_-dj_[iSequence]
2192                -newTolerance;
2193              //modification = max(modification,-dualTolerance_);
2194              //assert (fabs(modification)<1.0e-7);
2195              dj_[iSequence] += modification;
2196              cost_[iSequence] +=  modification;
2197              //cost_[iSequence+costOffset] += modification; // save change
2198#endif
2199            }
2200          }
2201        }
2202      }
2203    }
2204  }
2205
2206  if (sequenceIn_>=0) {
2207    lowerIn_ = lower_[sequenceIn_];
2208    upperIn_ = upper_[sequenceIn_];
2209    valueIn_ = solution_[sequenceIn_];
2210    dualIn_ = dj_[sequenceIn_];
2211
2212    if (numberTimesOptimal_) {
2213      // can we adjust cost back closer to original
2214      //*** add coding
2215    }
2216#if MODIFYCOST>1   
2217    // modify cost to hit zero exactly
2218    // so (dualIn_+modification)==theta_*alpha_
2219    double modification = theta_*alpha_-dualIn_;
2220    dualIn_ += modification;
2221    dj_[sequenceIn_]=dualIn_;
2222    cost_[sequenceIn_] += modification;
2223    //int costOffset = numberRows_+numberColumns_;
2224    //cost_[sequenceIn_+costOffset] += modification; // save change
2225    //assert (fabs(modification)<1.0e-6);
2226#ifdef CLP_DEBUG
2227    if ((handler_->logLevel()&32)&&fabs(modification)>1.0e-15)
2228      printf("exact %d new cost %g, change %g\n",sequenceIn_,
2229             cost_[sequenceIn_],modification);
2230#endif
2231#endif
2232   
2233    if (alpha_<0.0) {
2234      // as if from upper bound
2235      directionIn_=-1;
2236      upperIn_=valueIn_;
2237    } else {
2238      // as if from lower bound
2239      directionIn_=1;
2240      lowerIn_=valueIn_;
2241    }
2242  }
2243  //if (thisIncrease)
2244  //dualTolerance_+= 1.0e-6*dblParam_[ClpDualTolerance];
2245
2246  // clear arrays
2247
2248  for (i=0;i<2;i++) {
2249    memset(array[i],0,marker[i][0]*sizeof(double));
2250    memset(array[i]+marker[i][1],0,
2251           (numberColumns_-marker[i][1])*sizeof(double));
2252  }
2253}
2254/* Checks if tentative optimal actually means unbounded
2255   Returns -3 if not, 2 if is unbounded */
2256int 
2257ClpSimplexDual::checkUnbounded(CoinIndexedVector * ray,
2258                                   CoinIndexedVector * spare,
2259                                   double changeCost)
2260{
2261  int status=2; // say unbounded
2262  factorization_->updateColumn(spare,ray);
2263  // get reduced cost
2264  int i;
2265  int number=ray->getNumElements();
2266  int * index = ray->getIndices();
2267  double * array = ray->denseVector();
2268  for (i=0;i<number;i++) {
2269    int iRow=index[i];
2270    int iPivot=pivotVariable_[iRow];
2271    changeCost -= cost(iPivot)*array[iRow];
2272  }
2273  double way;
2274  if (changeCost>0.0) {
2275    //try going down
2276    way=1.0;
2277  } else if (changeCost<0.0) {
2278    //try going up
2279    way=-1.0;
2280  } else {
2281#ifdef CLP_DEBUG
2282    printf("can't decide on up or down\n");
2283#endif
2284    way=0.0;
2285    status=-3;
2286  }
2287  double movement=1.0e10*way; // some largish number
2288  double zeroTolerance = 1.0e-14*dualBound_;
2289  for (i=0;i<number;i++) {
2290    int iRow=index[i];
2291    int iPivot=pivotVariable_[iRow];
2292    double arrayValue = array[iRow];
2293    if (fabs(arrayValue)<zeroTolerance)
2294      arrayValue=0.0;
2295    double newValue=solution(iPivot)+movement*arrayValue;
2296    if (newValue>upper(iPivot)+primalTolerance_||
2297        newValue<lower(iPivot)-primalTolerance_)
2298      status=-3; // not unbounded
2299  }
2300  if (status==2) {
2301    // create ray
2302    delete [] ray_;
2303    ray_ = new double [numberColumns_];
2304    ClpFillN(ray_,numberColumns_,0.0);
2305    for (i=0;i<number;i++) {
2306      int iRow=index[i];
2307      int iPivot=pivotVariable_[iRow];
2308      double arrayValue = array[iRow];
2309      if (iPivot<numberColumns_&&fabs(arrayValue)>=zeroTolerance)
2310        ray_[iPivot] = way* array[iRow];
2311    }
2312  }
2313  ray->clear();
2314  return status;
2315}
2316/* Checks if finished.  Updates status */
2317void 
2318ClpSimplexDual::statusOfProblemInDual(int & lastCleaned,int type,
2319                                      ClpSimplexProgress * progress,
2320                                      double * givenDuals)
2321{
2322  bool normalType=true;
2323  if (type==2) {
2324    // trouble - restore solution
2325    memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
2326    memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
2327           numberRows_*sizeof(double));
2328    memcpy(columnActivityWork_,savedSolution_ ,
2329           numberColumns_*sizeof(double));
2330    forceFactorization_=1; // a bit drastic but ..
2331    changeMade_++; // say something changed
2332  }
2333  int tentativeStatus = problemStatus_;
2334  double changeCost;
2335  bool unflagVariables = true;
2336  if (problemStatus_>-3||factorization_->pivots()) {
2337    // factorize
2338    // later on we will need to recover from singularities
2339    // also we could skip if first time
2340    // save dual weights
2341    dualRowPivot_->saveWeights(this,1);
2342    if (type) {
2343      // is factorization okay?
2344      if (internalFactorize(1)) {
2345        // no - restore previous basis
2346        unflagVariables = false;
2347        assert (type==1);
2348        changeMade_++; // say something changed
2349        memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
2350        memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
2351               numberRows_*sizeof(double));
2352        memcpy(columnActivityWork_,savedSolution_ ,
2353               numberColumns_*sizeof(double));
2354        // get correct bounds on all variables
2355        double dummyChangeCost=0.0;
2356        changeBounds(true,rowArray_[2],dummyChangeCost);
2357        // throw away change
2358        int i;
2359        for (i=0;i<4;i++) 
2360          rowArray_[i]->clear();
2361        // need to reject something
2362        char x = isColumn(sequenceOut_) ? 'C' :'R';
2363        handler_->message(CLP_SIMPLEX_FLAG,messages_)
2364          <<x<<sequenceWithin(sequenceOut_)
2365          <<CoinMessageEol;
2366        setFlagged(sequenceOut_);
2367       
2368        forceFactorization_=1; // a bit drastic but ..
2369        type = 2;
2370        //assert (internalFactorize(1)==0);
2371        if (internalFactorize(1)) {
2372          memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
2373          memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
2374                 numberRows_*sizeof(double));
2375          memcpy(columnActivityWork_,savedSolution_ ,
2376                 numberColumns_*sizeof(double));
2377          // debug
2378          assert (internalFactorize(1)==0);
2379        }
2380      }
2381    }
2382    if (problemStatus_!=-4||factorization_->pivots()>10)
2383      problemStatus_=-3;
2384  }
2385  // at this stage status is -3 or -4 if looks infeasible
2386  // get primal and dual solutions
2387  gutsOfSolution(givenDuals,NULL);
2388  // Double check infeasibility if no action
2389  if (progress->lastIterationNumber(0)==numberIterations_) {
2390    if (dualRowPivot_->looksOptimal()) {
2391      numberPrimalInfeasibilities_ = 0;
2392      sumPrimalInfeasibilities_ = 0.0;
2393    }
2394  }
2395  // Check if looping
2396  int loop;
2397  if (!givenDuals&&type!=2) 
2398    loop = progress->looping();
2399  else
2400    loop=-1;
2401  int situationChanged=0;
2402  if (loop>=0) {
2403    problemStatus_ = loop; //exit if in loop
2404    if (!problemStatus_) {
2405      // declaring victory
2406      numberPrimalInfeasibilities_ = 0;
2407      sumPrimalInfeasibilities_ = 0.0;
2408    }
2409    return;
2410  } else if (loop<-1) {
2411    // something may have changed
2412    gutsOfSolution(NULL,NULL);
2413    situationChanged=1;
2414  }
2415  // really for free variables in
2416  if((progressFlag_&2)!=0) {
2417    situationChanged=2;
2418  }
2419  progressFlag_ = 0; //reset progress flag
2420  handler_->message(CLP_SIMPLEX_STATUS,messages_)
2421    <<numberIterations_<<objectiveValue();
2422  handler_->printing(sumPrimalInfeasibilities_>0.0)
2423    <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
2424  handler_->printing(sumDualInfeasibilities_>0.0)
2425    <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
2426  handler_->printing(numberDualInfeasibilitiesWithoutFree_
2427                     <numberDualInfeasibilities_)
2428                       <<numberDualInfeasibilitiesWithoutFree_;
2429  handler_->message()<<CoinMessageEol;
2430
2431  // dual bound coming in
2432  //double saveDualBound = dualBound_;
2433  while (problemStatus_<=-3) {
2434    int cleanDuals=0;
2435    if (situationChanged!=0)
2436      cleanDuals=1;
2437    int numberChangedBounds=0;
2438    int doOriginalTolerance=0;
2439    if ( lastCleaned==numberIterations_)
2440      doOriginalTolerance=1;
2441    // check optimal
2442    // give code benefit of doubt
2443    if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
2444        sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
2445      // say optimal (with these bounds etc)
2446      numberDualInfeasibilities_ = 0;
2447      sumDualInfeasibilities_ = 0.0;
2448      numberPrimalInfeasibilities_ = 0;
2449      sumPrimalInfeasibilities_ = 0.0;
2450    }
2451    //if (dualFeasible()||problemStatus_==-4||(primalFeasible()&&!numberDualInfeasibilitiesWithoutFree_)) {
2452    if (dualFeasible()||problemStatus_==-4) {
2453      progress->modifyObjective(objectiveValue_
2454                               -sumDualInfeasibilities_*dualBound_);
2455      if (primalFeasible()) {
2456        normalType=false;
2457        // may be optimal - or may be bounds are wrong
2458        handler_->message(CLP_DUAL_BOUNDS,messages_)
2459          <<dualBound_
2460          <<CoinMessageEol;
2461        // save solution in case unbounded
2462        ClpDisjointCopyN(columnActivityWork_,numberColumns_,
2463                          columnArray_[0]->denseVector());
2464        ClpDisjointCopyN(rowActivityWork_,numberRows_,
2465                          rowArray_[2]->denseVector());
2466        numberChangedBounds=changeBounds(false,rowArray_[3],changeCost);
2467        if (numberChangedBounds<=0&&!numberDualInfeasibilities_) {
2468          //looks optimal - do we need to reset tolerance
2469          if (perturbation_==101) {
2470            perturbation_=102; // stop any perturbations
2471            cleanDuals=1;
2472            createRim(4);
2473          }
2474          if (lastCleaned<numberIterations_&&numberTimesOptimal_<4) {
2475            doOriginalTolerance=2;
2476            numberTimesOptimal_++;
2477            changeMade_++; // say something changed
2478            if (numberTimesOptimal_==1) {
2479              dualTolerance_ = dblParam_[ClpDualTolerance];
2480              // better to have small tolerance even if slower
2481              factorization_->zeroTolerance(1.0e-15);
2482            } else {
2483              dualTolerance_ = dblParam_[ClpDualTolerance];
2484              dualTolerance_ *= pow(2.0,numberTimesOptimal_-1);
2485            }
2486            cleanDuals=2; // If nothing changed optimal else primal
2487          } else {
2488            problemStatus_=0; // optimal
2489            if (lastCleaned<numberIterations_) {
2490              handler_->message(CLP_SIMPLEX_GIVINGUP,messages_)
2491                <<CoinMessageEol;
2492            }
2493          }
2494        } else {
2495          cleanDuals=1;
2496          if (doOriginalTolerance==1) {
2497            // check unbounded
2498            // find a variable with bad dj
2499            int iSequence;
2500            int iChosen=-1;
2501            double largest = 100.0*primalTolerance_;
2502            for (iSequence=0;iSequence<numberRows_+numberColumns_;
2503                 iSequence++) {
2504              double djValue = dj_[iSequence];
2505              double originalLo = originalLower(iSequence);
2506              double originalUp = originalUpper(iSequence);
2507              if (fabs(djValue)>fabs(largest)) {
2508                if (getStatus(iSequence)!=basic) {
2509                  if (djValue>0&&originalLo<-1.0e20) {
2510                    if (djValue>fabs(largest)) {
2511                      largest=djValue;
2512                      iChosen=iSequence;
2513                    }
2514                  } else if (djValue<0&&originalUp>1.0e20) {
2515                    if (-djValue>fabs(largest)) {
2516                      largest=djValue;
2517                      iChosen=iSequence;
2518                    }
2519                  } 
2520                }
2521              }
2522            }
2523            if (iChosen>=0) {
2524              int iSave=sequenceIn_;
2525              sequenceIn_=iChosen;
2526              unpack(rowArray_[1]);
2527              sequenceIn_ = iSave;
2528              // if dual infeasibilities then must be free vector so add in dual
2529              if (numberDualInfeasibilities_) {
2530                if (fabs(changeCost)>1.0e-5)
2531                  printf("Odd free/unbounded combo\n");
2532                changeCost += cost_[iChosen];
2533              }
2534              problemStatus_ = checkUnbounded(rowArray_[1],rowArray_[0],
2535                                              changeCost);
2536              rowArray_[1]->clear();
2537            } else {
2538              problemStatus_=-3;
2539            }
2540            if (problemStatus_==2&&perturbation_==101) {
2541              perturbation_=102; // stop any perturbations
2542              cleanDuals=1;
2543              createRim(4);
2544              problemStatus_=-1;
2545            }
2546            if (problemStatus_==2) {
2547              // it is unbounded - restore solution
2548              // but first add in changes to non-basic
2549              int iColumn;
2550              double * original = columnArray_[0]->denseVector();
2551              for (iColumn=0;iColumn<numberColumns_;iColumn++) {
2552                if(getColumnStatus(iColumn)!= basic)
2553                  ray_[iColumn] += 
2554                    columnActivityWork_[iColumn]-original[iColumn];
2555                columnActivityWork_[iColumn] = original[iColumn];
2556              }
2557              ClpDisjointCopyN(rowArray_[2]->denseVector(),numberRows_,
2558                                rowActivityWork_);
2559            }
2560          } else {
2561            doOriginalTolerance=2;
2562            rowArray_[0]->clear();
2563          }
2564        }
2565        ClpFillN(columnArray_[0]->denseVector(),numberColumns_,0.0);
2566        ClpFillN(rowArray_[2]->denseVector(),numberRows_,0.0);
2567      } 
2568      if (problemStatus_==-4) {
2569        // may be infeasible - or may be bounds are wrong
2570        numberChangedBounds=changeBounds(false,NULL,changeCost);
2571        if (perturbation_==101) {
2572          perturbation_=102; // stop any perturbations
2573          cleanDuals=1;
2574          numberChangedBounds=1;
2575          createRim(4);
2576        }
2577        if (numberChangedBounds<=0||dualBound_>1.0e20||
2578            (largestPrimalError_>1.0&&dualBound_>1.0e17)) {
2579          problemStatus_=1; // infeasible
2580        } else {
2581          normalType=false;
2582          problemStatus_=-1; //iterate
2583          cleanDuals=1;
2584          if (numberChangedBounds<=0)
2585            doOriginalTolerance=2;
2586          // and delete ray which has been created
2587          delete [] ray_;
2588          ray_ = NULL;
2589        }
2590
2591      }
2592    } else {
2593      cleanDuals=1;
2594    }
2595    if (problemStatus_<0) {
2596      if (doOriginalTolerance==2) {
2597        // put back original tolerance
2598        lastCleaned=numberIterations_;
2599        handler_->message(CLP_DUAL_ORIGINAL,messages_)
2600          <<CoinMessageEol;
2601
2602        perturbation_=102; // stop any perturbations
2603#if 0
2604        double * xcost = new double[numberRows_+numberColumns_];
2605        double * xlower = new double[numberRows_+numberColumns_];
2606        double * xupper = new double[numberRows_+numberColumns_];
2607        double * xdj = new double[numberRows_+numberColumns_];
2608        double * xsolution = new double[numberRows_+numberColumns_];
2609        memcpy(xcost,cost_,(numberRows_+numberColumns_)*sizeof(double));
2610        memcpy(xlower,lower_,(numberRows_+numberColumns_)*sizeof(double));
2611        memcpy(xupper,upper_,(numberRows_+numberColumns_)*sizeof(double));
2612        memcpy(xdj,dj_,(numberRows_+numberColumns_)*sizeof(double));
2613        memcpy(xsolution,solution_,(numberRows_+numberColumns_)*sizeof(double));
2614#endif
2615        createRim(4);
2616        // make sure duals are current
2617        computeDuals(givenDuals);
2618        checkDualSolution();
2619#if 0
2620        int i;
2621        for (i=0;i<numberRows_+numberColumns_;i++) {
2622          if (cost_[i]!=xcost[i])
2623            printf("** %d old cost %g new %g sol %g\n",
2624                   i,xcost[i],cost_[i],solution_[i]);
2625          if (lower_[i]!=xlower[i])
2626            printf("** %d old lower %g new %g sol %g\n",
2627                   i,xlower[i],lower_[i],solution_[i]);
2628          if (upper_[i]!=xupper[i])
2629            printf("** %d old upper %g new %g sol %g\n",
2630                   i,xupper[i],upper_[i],solution_[i]);
2631          if (dj_[i]!=xdj[i])
2632            printf("** %d old dj %g new %g sol %g\n",
2633                   i,xdj[i],dj_[i],solution_[i]);
2634          if (solution_[i]!=xsolution[i])
2635            printf("** %d old solution %g new %g sol %g\n",
2636                   i,xsolution[i],solution_[i],solution_[i]);
2637        }
2638        //delete [] xcost;
2639        //delete [] xupper;
2640        //delete [] xlower;
2641        //delete [] xdj;
2642        //delete [] xsolution;
2643#endif
2644        // put back bounds as they were if was optimal
2645        if (doOriginalTolerance==2) {
2646          changeMade_++; // say something changed
2647          changeBounds(true,NULL,changeCost);
2648          cleanDuals=2;
2649          //cleanDuals=1;
2650        }
2651#if 0
2652        //int i;
2653        for (i=0;i<numberRows_+numberColumns_;i++) {
2654          if (cost_[i]!=xcost[i])
2655            printf("** %d old cost %g new %g sol %g\n",
2656                   i,xcost[i],cost_[i],solution_[i]);
2657          if (lower_[i]!=xlower[i])
2658            printf("** %d old lower %g new %g sol %g\n",
2659                   i,xlower[i],lower_[i],solution_[i]);
2660          if (upper_[i]!=xupper[i])
2661            printf("** %d old upper %g new %g sol %g\n",
2662                   i,xupper[i],upper_[i],solution_[i]);
2663          if (dj_[i]!=xdj[i])
2664            printf("** %d old dj %g new %g sol %g\n",
2665                   i,xdj[i],dj_[i],solution_[i]);
2666          if (solution_[i]!=xsolution[i])
2667            printf("** %d old solution %g new %g sol %g\n",
2668                   i,xsolution[i],solution_[i],solution_[i]);
2669        }
2670        delete [] xcost;
2671        delete [] xupper;
2672        delete [] xlower;
2673        delete [] xdj;
2674        delete [] xsolution;
2675#endif
2676      }
2677      if (cleanDuals==1||(cleanDuals==2&&!numberDualInfeasibilities_)) {
2678        // make sure dual feasible
2679        // look at all rows and columns
2680        rowArray_[0]->clear();
2681        columnArray_[0]->clear();
2682        double objectiveChange=0.0;
2683#if 0
2684        double * xcost = new double[numberRows_+numberColumns_];
2685        double * xlower = new double[numberRows_+numberColumns_];
2686        double * xupper = new double[numberRows_+numberColumns_];
2687        double * xdj = new double[numberRows_+numberColumns_];
2688        double * xsolution = new double[numberRows_+numberColumns_];
2689        memcpy(xcost,cost_,(numberRows_+numberColumns_)*sizeof(double));
2690        memcpy(xlower,lower_,(numberRows_+numberColumns_)*sizeof(double));
2691        memcpy(xupper,upper_,(numberRows_+numberColumns_)*sizeof(double));
2692        memcpy(xdj,dj_,(numberRows_+numberColumns_)*sizeof(double));
2693        memcpy(xsolution,solution_,(numberRows_+numberColumns_)*sizeof(double));
2694#endif
2695        updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[1],
2696                          0.0,objectiveChange,true);
2697#if 0
2698        int i;
2699        for (i=0;i<numberRows_+numberColumns_;i++) {
2700          if (cost_[i]!=xcost[i])
2701            printf("** %d old cost %g new %g sol %g\n",
2702                   i,xcost[i],cost_[i],solution_[i]);
2703          if (lower_[i]!=xlower[i])
2704            printf("** %d old lower %g new %g sol %g\n",
2705                   i,xlower[i],lower_[i],solution_[i]);
2706          if (upper_[i]!=xupper[i])
2707            printf("** %d old upper %g new %g sol %g\n",
2708                   i,xupper[i],upper_[i],solution_[i]);
2709          if (dj_[i]!=xdj[i])
2710            printf("** %d old dj %g new %g sol %g\n",
2711                   i,xdj[i],dj_[i],solution_[i]);
2712          if (solution_[i]!=xsolution[i])
2713            printf("** %d old solution %g new %g sol %g\n",
2714                   i,xsolution[i],solution_[i],solution_[i]);
2715        }
2716        delete [] xcost;
2717        delete [] xupper;
2718        delete [] xlower;
2719        delete [] xdj;
2720        delete [] xsolution;
2721#endif
2722        // for now - recompute all
2723        gutsOfSolution(NULL,NULL);
2724        //assert(numberDualInfeasibilitiesWithoutFree_==0);
2725
2726        if (numberDualInfeasibilities_||situationChanged==2) 
2727          problemStatus_=-1; // carry on as normal
2728        situationChanged=0;
2729      } else {
2730        // iterate
2731        if (cleanDuals!=2) 
2732          problemStatus_=-1;
2733        else
2734          problemStatus_ = 10; // try primal
2735      }
2736    }
2737  }
2738  if (type==0||type==1) {
2739    if (!type) {
2740      // create save arrays
2741      delete [] saveStatus_;
2742      delete [] savedSolution_;
2743      saveStatus_ = new unsigned char [numberRows_+numberColumns_];
2744      savedSolution_ = new double [numberRows_+numberColumns_];
2745    }
2746    // save arrays
2747    memcpy(saveStatus_,status_,(numberColumns_+numberRows_)*sizeof(char));
2748    memcpy(savedSolution_+numberColumns_ ,rowActivityWork_,
2749           numberRows_*sizeof(double));
2750    memcpy(savedSolution_ ,columnActivityWork_,numberColumns_*sizeof(double));
2751  }
2752
2753  // restore weights (if saved) - also recompute infeasibility list
2754  if (tentativeStatus>-3) 
2755    dualRowPivot_->saveWeights(this,(type <2) ? 2 : 4);
2756  else
2757    dualRowPivot_->saveWeights(this,3);
2758  // unflag all variables (we may want to wait a bit?)
2759  if (tentativeStatus!=-2&&unflagVariables) {
2760    int iRow;
2761    for (iRow=0;iRow<numberRows_;iRow++) {
2762      int iPivot=pivotVariable_[iRow];
2763      clearFlagged(iPivot);
2764    }
2765  }
2766  // see if cutoff reached
2767  double limit = 0.0;
2768  getDblParam(ClpDualObjectiveLimit, limit);
2769  if(fabs(limit)<1.0e30&&objectiveValue()*optimizationDirection_>
2770           optimizationDirection_*limit&&
2771           !numberAtFakeBound()&&!numberDualInfeasibilities_) {
2772    problemStatus_=1;
2773    secondaryStatus_ = 1; // and say was on cutoff
2774  }
2775  if (problemStatus_<0&&!changeMade_) {
2776    problemStatus_=4; // unknown
2777  }
2778  lastGoodIteration_ = numberIterations_;
2779#if 0
2780  double thisObj = progress->lastObjective(0);
2781  double lastObj = progress->lastObjective(1);
2782  if (lastObj>thisObj+1.0e-6*max(fabs(thisObj),fabs(lastObj))+1.0e-8
2783      &&givenDuals==NULL) {
2784    int maxFactor = factorization_->maximumPivots();
2785    if (maxFactor>10) {
2786      if (forceFactorization_<0)
2787        forceFactorization_= maxFactor;
2788      forceFactorization_ = max (1,(forceFactorization_>>1));
2789      printf("Reducing factorization frequency\n");
2790    }
2791  }
2792#endif
2793}
2794/* While updateDualsInDual sees what effect is of flip
2795   this does actual flipping.
2796   If change >0.0 then value in array >0.0 => from lower to upper
2797*/
2798void 
2799ClpSimplexDual::flipBounds(CoinIndexedVector * rowArray,
2800                  CoinIndexedVector * columnArray,
2801                  double change)
2802{
2803  int number;
2804  int * which;
2805 
2806  int iSection;
2807
2808  for (iSection=0;iSection<2;iSection++) {
2809    int i;
2810    double * solution = solutionRegion(iSection);
2811    double * lower = lowerRegion(iSection);
2812    double * upper = upperRegion(iSection);
2813    int addSequence;
2814    if (!iSection) {
2815      number = rowArray->getNumElements();
2816      which = rowArray->getIndices();
2817      addSequence = numberColumns_;
2818    } else {
2819      number = columnArray->getNumElements();
2820      which = columnArray->getIndices();
2821      addSequence = 0;
2822    }
2823   
2824    for (i=0;i<number;i++) {
2825      int iSequence = which[i];
2826      Status status = getStatus(iSequence+addSequence);
2827
2828      switch(status) {
2829
2830      case basic:
2831      case isFree:
2832      case superBasic:
2833      case ClpSimplex::isFixed:
2834        break;
2835      case atUpperBound:
2836        // to lower bound
2837        setStatus(iSequence+addSequence,atLowerBound);
2838        solution[iSequence] = lower[iSequence];
2839        break;
2840      case atLowerBound:
2841        // to upper bound
2842        setStatus(iSequence+addSequence,atUpperBound);
2843        solution[iSequence] = upper[iSequence];
2844        break;
2845      }
2846    }
2847  }
2848  rowArray->setNumElements(0);
2849  columnArray->setNumElements(0);
2850}
2851// Restores bound to original bound
2852void 
2853ClpSimplexDual::originalBound( int iSequence)
2854{
2855  if (getFakeBound(iSequence)!=noFake)
2856    numberFake_--;;
2857  if (iSequence>=numberColumns_) {
2858    // rows
2859    int iRow = iSequence-numberColumns_;
2860    rowLowerWork_[iRow]=rowLower_[iRow];
2861    rowUpperWork_[iRow]=rowUpper_[iRow];
2862    if (rowScale_) {
2863      if (rowLowerWork_[iRow]>-1.0e50)
2864        rowLowerWork_[iRow] *= rowScale_[iRow];
2865      if (rowUpperWork_[iRow]<1.0e50)
2866        rowUpperWork_[iRow] *= rowScale_[iRow];
2867    }
2868  } else {
2869    // columns
2870    columnLowerWork_[iSequence]=columnLower_[iSequence];
2871    columnUpperWork_[iSequence]=columnUpper_[iSequence];
2872    if (rowScale_) {
2873      double multiplier = 1.0/columnScale_[iSequence];
2874      if (columnLowerWork_[iSequence]>-1.0e50)
2875        columnLowerWork_[iSequence] *= multiplier;
2876      if (columnUpperWork_[iSequence]<1.0e50)
2877        columnUpperWork_[iSequence] *= multiplier;
2878    }
2879  }
2880  setFakeBound(iSequence,noFake);
2881}
2882/* As changeBounds but just changes new bounds for a single variable.
2883   Returns true if change */
2884bool 
2885ClpSimplexDual::changeBound( int iSequence)
2886{
2887  // old values
2888  double oldLower=lower_[iSequence];
2889  double oldUpper=upper_[iSequence];
2890  double value=solution_[iSequence];
2891  bool modified=false;
2892  originalBound(iSequence);
2893  // original values
2894  double lowerValue=lower_[iSequence];
2895  double upperValue=upper_[iSequence];
2896  // back to altered values
2897  lower_[iSequence] = oldLower;
2898  upper_[iSequence] = oldUpper;
2899  if (getFakeBound(iSequence)!=noFake)
2900    numberFake_--;;
2901  if (value==oldLower) {
2902    if (upperValue > oldLower + dualBound_) {
2903      upper_[iSequence]=oldLower+dualBound_;
2904      setFakeBound(iSequence,upperFake);
2905      modified=true;
2906      numberFake_++;
2907    }
2908  } else if (value==oldUpper) {
2909    if (lowerValue < oldUpper - dualBound_) {
2910      lower_[iSequence]=oldUpper-dualBound_;
2911      setFakeBound(iSequence,lowerFake);
2912      modified=true;
2913      numberFake_++;
2914    }
2915  } else {
2916    assert(value==oldLower||value==oldUpper);
2917  }
2918  return modified;
2919}
2920// Perturbs problem
2921void 
2922ClpSimplexDual::perturb()
2923{
2924  if (perturbation_>100)
2925    return; //perturbed already
2926  if (perturbation_==100)
2927    perturbation_=50; // treat as normal
2928  bool modifyRowCosts=false;
2929  // dual perturbation
2930  double perturbation=1.0e-20;
2931  // maximum fraction of cost to perturb
2932  double maximumFraction = 1.0e-5;
2933  double constantPerturbation = 100.0*dualTolerance_;
2934  const int * lengths = matrix_->getVectorLengths();
2935  int maxLength=0;
2936  int minLength=numberRows_;
2937  if (!numberIterations_&&perturbation_==50) {
2938    // See if we need to perturb
2939    double * sort = new double[numberColumns_];
2940    // Use objective BEFORE scaling
2941    const double * obj = objective();
2942    int i;
2943    for (i=0;i<numberColumns_;i++) {
2944      double value = fabs(obj[i]);
2945      sort[i]=value;
2946    }
2947    std::sort(sort,sort+numberColumns_);
2948    int number=1;
2949    double last = sort[0];
2950    for (i=1;i<numberColumns_;i++) {
2951      if (last!=sort[i])
2952        number++;
2953      last=sort[i];
2954    }
2955    delete [] sort;
2956#if 0
2957    printf("nnz %d percent %d",number,(number*100)/numberColumns_);
2958    if (number*4>numberColumns_)
2959      printf(" - Would not perturb\n");
2960    else
2961      printf(" - Would perturb\n");
2962    exit(0);
2963#endif
2964    if (number*4>numberColumns_) {
2965      perturbation_=100;
2966      return; // good enough
2967    }
2968  }
2969  int iColumn;
2970  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
2971    if (columnLowerWork_[iColumn]<columnUpperWork_[iColumn]) {
2972      int length = lengths[iColumn];
2973      if (length>2) {
2974        maxLength = max(maxLength,length);
2975        minLength = min(minLength,length);
2976      }
2977    }
2978  }
2979  // If > 70 then do rows
2980  if (perturbation_>=70) {
2981    modifyRowCosts=true;
2982    perturbation_ -= 20;
2983    printf("Row costs modified, ");
2984  }
2985  bool uniformChange=false;
2986  if (perturbation_>50) {
2987    // Experiment
2988    // maximumFraction could be 1.0e-10 to 1.0
2989    double m[]={1.0e-10,1.0e-9,1.0e-8,1.0e-7,1.0e-6,1.0e-5,1.0e-4,1.0e-3,1.0e-2,1.0e-1,1.0};
2990    maximumFraction = m[min(perturbation_-51,10)];
2991  }
2992  int iRow;
2993  double smallestNonZero=1.0e100;
2994  int numberNonZero=0;
2995  if (perturbation_>=50) {
2996    perturbation = 1.0e-8;
2997    bool allSame=true;
2998    double lastValue=0.0;
2999    for (iRow=0;iRow<numberRows_;iRow++) {
3000      double lo = rowLowerWork_[iRow];
3001      double up = rowUpperWork_[iRow];
3002      if (lo<up) {
3003        double value = fabs(rowObjectiveWork_[iRow]);
3004        perturbation = max(perturbation,value);
3005        if (value) {
3006          modifyRowCosts=true;
3007          smallestNonZero = min(smallestNonZero,value);
3008        }
3009      } 
3010      if (lo&&lo>-1.0e10) {
3011        numberNonZero++;
3012        lo=fabs(lo);
3013        if (!lastValue) 
3014          lastValue=lo;
3015        else if (fabs(lo-lastValue)>1.0e-7)
3016          allSame=false;
3017      }
3018      if (up&&up<1.0e10) {
3019        numberNonZero++;
3020        up=fabs(up);
3021        if (!lastValue) 
3022          lastValue=up;
3023        else if (fabs(up-lastValue)>1.0e-7)
3024          allSame=false;
3025      }
3026    }
3027    double lastValue2=0.0;
3028    for (iColumn=0;iColumn<numberColumns_;iColumn++) { 
3029      double lo = columnLowerWork_[iColumn];
3030      double up = columnUpperWork_[iColumn];
3031      if (lo<up) {
3032        double value = 
3033          fabs(objectiveWork_[iColumn]);
3034        perturbation = max(perturbation,value);
3035        if (value) {
3036          smallestNonZero = min(smallestNonZero,value);
3037        }
3038      }
3039      if (lo&&lo>-1.0e10) {
3040        //numberNonZero++;
3041        lo=fabs(lo);
3042        if (!lastValue2) 
3043          lastValue2=lo;
3044        else if (fabs(lo-lastValue2)>1.0e-7)
3045          allSame=false;
3046      }
3047      if (up&&up<1.0e10) {
3048        //numberNonZero++;
3049        up=fabs(up);
3050        if (!lastValue2) 
3051          lastValue2=up;
3052        else if (fabs(up-lastValue2)>1.0e-7)
3053          allSame=false;
3054      }
3055    }
3056    if (allSame) {
3057      // Check elements
3058      double smallestNegative;
3059      double largestNegative;
3060      double smallestPositive;
3061      double largestPositive;
3062      matrix_->rangeOfElements(smallestNegative,largestNegative,
3063                               smallestPositive,largestPositive);
3064      if (smallestNegative==largestNegative&&
3065          smallestPositive==largestPositive) {
3066        // Really hit perturbation
3067        maximumFraction = max(1.0e-3*max(lastValue,lastValue2),maximumFraction);
3068      }
3069    }
3070    perturbation = min(perturbation,smallestNonZero/maximumFraction);
3071  } else {
3072    // user is in charge
3073    maximumFraction = 1.0e-1;
3074    // but some experiments
3075    if (perturbation_<=-900) {
3076      modifyRowCosts=true;
3077      perturbation_ += 1000;
3078      printf("Row costs modified, ");
3079    }
3080    if (perturbation_<=-10) {
3081      perturbation_ += 10; 
3082      maximumFraction = 1.0;
3083      if ((-perturbation_)%100>=10) {
3084        uniformChange=true;
3085        perturbation_+=20;
3086      }
3087      while (perturbation_<-10) {
3088        perturbation_ += 100;
3089        maximumFraction *= 1.0e-1;
3090      }
3091    }
3092    perturbation = pow(10.0,perturbation_);
3093  }
3094  double largestZero=0.0;
3095  double largest=0.0;
3096  double largestPerCent=0.0;
3097  // modify costs
3098  bool printOut=(handler_->logLevel()==63);
3099  printOut=false;
3100  if (modifyRowCosts) {
3101    for (iRow=0;iRow<numberRows_;iRow++) {
3102      if (rowLowerWork_[iRow]<rowUpperWork_[iRow]) {
3103        double value = perturbation;
3104        double currentValue = rowObjectiveWork_[iRow];
3105        value = min(value,maximumFraction*(fabs(currentValue)+1.0e-1*perturbation+1.0e-3));
3106        if (rowLowerWork_[iRow]>-largeValue_) {
3107          if (fabs(rowLowerWork_[iRow])<fabs(rowUpperWork_[iRow])) 
3108            value *= CoinDrand48();
3109          else
3110            value *= -CoinDrand48();
3111        } else if (rowUpperWork_[iRow]<largeValue_) {
3112          value *= -CoinDrand48();
3113        } else {
3114          value=0.0;
3115        }
3116        if (currentValue) {
3117          largest = max(largest,fabs(value));
3118          if (fabs(value)>fabs(currentValue)*largestPerCent)
3119            largestPerCent=fabs(value/currentValue);
3120        } else {
3121          largestZero=max(largestZero,fabs(value));
3122        }
3123        if (printOut)
3124          printf("row %d cost %g change %g\n",iRow,rowObjectiveWork_[iRow],value);
3125        rowObjectiveWork_[iRow] += value;
3126      }
3127    }
3128  }
3129  // more its but faster double weight[]={1.0e-4,1.0e-2,1.0e-1,1.0,2.0,10.0,100.0,200.0,400.0,600.0,1000.0};
3130  // good its double weight[]={1.0e-4,1.0e-2,5.0e-1,1.0,2.0,5.0,10.0,20.0,30.0,40.0,100.0};
3131  double weight[]={1.0e-4,1.0e-2,5.0e-1,1.0,2.0,5.0,10.0,20.0,30.0,40.0,100.0};
3132  //double weight[]={1.0e-4,1.0e-2,5.0e-1,1.0,20.0,50.0,100.0,120.0,130.0,140.0,200.0};
3133  double extraWeight=10.0;
3134  // Scale back if wanted
3135  double weight2[]={1.0e-4,1.0e-2,5.0e-1,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0};
3136  if (constantPerturbation<99.0*dualTolerance_) {
3137    perturbation *= 0.1;
3138    extraWeight=0.5;
3139    memcpy(weight,weight2,sizeof(weight2));
3140  }
3141  // adjust weights if all columns long
3142  double factor=1.0;
3143  if (maxLength) {
3144    factor = 3.0/(double) minLength;
3145  }
3146  // Make variables with more elements more expensive
3147  const double m1 = 0.5;
3148  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
3149    if (columnLowerWork_[iColumn]<columnUpperWork_[iColumn]&&getStatus(iColumn)!=basic) {
3150      double value = perturbation;
3151      double currentValue = objectiveWork_[iColumn];
3152      value = min(value,constantPerturbation+maximumFraction*(fabs(currentValue)+1.0e-1*perturbation+1.0e-8));
3153      //value = min(value,constantPerturbation;+maximumFraction*fabs(currentValue));
3154      double value2 = constantPerturbation+1.0e-1*smallestNonZero;
3155      if (uniformChange) {
3156        value = maximumFraction;
3157        value2=maximumFraction;
3158      }
3159      if (columnLowerWork_[iColumn]>-largeValue_) {
3160        if (fabs(columnLowerWork_[iColumn])<
3161            fabs(columnUpperWork_[iColumn])) {
3162          value *= (1.0-m1 +m1*CoinDrand48());
3163          value2 *= (1.0-m1 +m1*CoinDrand48());
3164        } else {
3165          value *= -(1.0-m1+m1*CoinDrand48());
3166          value2 *= -(1.0-m1+m1*CoinDrand48());
3167        }
3168      } else if (columnUpperWork_[iColumn]<largeValue_) {
3169        value *= -(1.0-m1+m1*CoinDrand48());
3170        value2 *= -(1.0-m1+m1*CoinDrand48());
3171      } else {
3172        value=0.0;
3173      }
3174      if (value) {
3175        int length = lengths[iColumn];
3176        if (length>3) {
3177          length = (int) ((double) length * factor);
3178          length = max(3,length);
3179        }
3180        double multiplier;
3181#if 1
3182        if (length<10)
3183          multiplier=weight[length];
3184        else
3185          multiplier = weight[10];
3186#else
3187        if (length<10)
3188          multiplier=weight[length];
3189        else
3190          multiplier = weight[10]+extraWeight*(length-10);
3191        multiplier *= 0.5;
3192#endif
3193        value *= multiplier;
3194        value = min (value,value2);
3195        if (fabs(value)<=dualTolerance_)
3196          value=0.0;
3197        if (currentValue) {
3198          largest = max(largest,fabs(value));
3199          if (fabs(value)>fabs(currentValue)*largestPerCent)
3200            largestPerCent=fabs(value/currentValue);
3201        } else {
3202          largestZero=max(largestZero,fabs(value));
3203        }
3204        if (printOut)
3205          printf("col %d cost %g change %g\n",iColumn,objectiveWork_[iColumn],value);
3206        objectiveWork_[iColumn] += value;
3207      }
3208    }
3209  }
3210  handler_->message(CLP_SIMPLEX_PERTURB,messages_)
3211    <<100.0*maximumFraction<<perturbation<<largest<<100.0*largestPerCent<<largestZero
3212    <<CoinMessageEol;
3213  // and zero changes
3214  //int nTotal = numberRows_+numberColumns_;
3215  //memset(cost_+nTotal,0,nTotal*sizeof(double));
3216  // say perturbed
3217  perturbation_=101;
3218
3219}
3220/* For strong branching.  On input lower and upper are new bounds
3221   while on output they are change in objective function values
3222   (>1.0e50 infeasible).
3223   Return code is 0 if nothing interesting, -1 if infeasible both
3224   ways and +1 if infeasible one way (check values to see which one(s))
3225*/
3226int ClpSimplexDual::strongBranching(int numberVariables,const int * variables,
3227                                    double * newLower, double * newUpper,
3228                                    double ** outputSolution,
3229                                    int * outputStatus, int * outputIterations,
3230                                    bool stopOnFirstInfeasible,
3231                                    bool alwaysFinish)
3232{
3233  int i;
3234  int returnCode=0;
3235  double saveObjectiveValue = objectiveValue_;
3236#if 1
3237  algorithm_ = -1;
3238
3239  //scaling(false);
3240
3241  // put in standard form (and make row copy)
3242  // create modifiable copies of model rim and do optional scaling
3243  createRim(7+8+16,true);
3244
3245  // change newLower and newUpper if scaled
3246
3247  // Do initial factorization
3248  // and set certain stuff
3249  // We can either set increasing rows so ...IsBasic gives pivot row
3250  // or we can just increment iBasic one by one
3251  // for now let ...iBasic give pivot row
3252  factorization_->increasingRows(2);
3253  // row activities have negative sign
3254  factorization_->slackValue(-1.0);
3255  factorization_->zeroTolerance(1.0e-13);
3256
3257  int factorizationStatus = internalFactorize(0);
3258  if (factorizationStatus<0)
3259    return 1; // some error
3260  else if (factorizationStatus)
3261    handler_->message(CLP_SINGULARITIES,messages_)
3262      <<factorizationStatus
3263      <<CoinMessageEol;
3264 
3265  // save stuff
3266  ClpFactorization saveFactorization(*factorization_);
3267  // save basis and solution
3268  double * saveSolution = new double[numberRows_+numberColumns_];
3269  memcpy(saveSolution,solution_,
3270         (numberRows_+numberColumns_)*sizeof(double));
3271  unsigned char * saveStatus = 
3272    new unsigned char [numberRows_+numberColumns_];
3273  memcpy(saveStatus,status_,(numberColumns_+numberRows_)*sizeof(char));
3274  // save bounds as createRim makes clean copies
3275  double * saveLower = new double[numberRows_+numberColumns_];
3276  memcpy(saveLower,lower_,
3277         (numberRows_+numberColumns_)*sizeof(double));
3278  double * saveUpper = new double[numberRows_+numberColumns_];
3279  memcpy(saveUpper,upper_,
3280         (numberRows_+numberColumns_)*sizeof(double));
3281  double * saveObjective = new double[numberRows_+numberColumns_];
3282  memcpy(saveObjective,cost_,
3283         (numberRows_+numberColumns_)*sizeof(double));
3284  int * savePivot = new int [numberRows_];
3285  memcpy(savePivot, pivotVariable_, numberRows_*sizeof(int));
3286  // need to save/restore weights.
3287
3288  int iSolution = 0;
3289  for (i=0;i<numberVariables;i++) {
3290    int iColumn = variables[i];
3291    double objectiveChange;
3292    double saveBound;
3293   
3294    // try down
3295   
3296    saveBound = columnUpper_[iColumn];
3297    // external view - in case really getting optimal
3298    columnUpper_[iColumn] = newUpper[i];
3299    if (scalingFlag_<=0) 
3300      upper_[iColumn] = newUpper[i];
3301    else 
3302      upper_[iColumn] = newUpper[i]/columnScale_[iColumn]; // scale
3303    // Start of fast iterations
3304    int status = fastDual(alwaysFinish);
3305    if (status) {
3306      // not finished - might be optimal
3307      checkPrimalSolution(rowActivityWork_,columnActivityWork_);
3308      double limit = 0.0;
3309      getDblParam(ClpDualObjectiveLimit, limit);
3310      if (!numberPrimalInfeasibilities_&&objectiveValue()*optimizationDirection_<
3311          optimizationDirection_*limit) { 
3312        problemStatus_=0;
3313      } 
3314      status=problemStatus_;
3315    }
3316
3317    if (scalingFlag_<=0) {
3318      memcpy(outputSolution[iSolution],solution_,numberColumns_*sizeof(double));
3319    } else {
3320      int j;
3321      double * sol = outputSolution[iSolution];
3322      for (j=0;j<numberColumns_;j++) 
3323        sol[j] = solution_[j]*columnScale_[j];
3324    }
3325    outputStatus[iSolution]=status;
3326    outputIterations[iSolution]=numberIterations_;
3327    iSolution++;
3328    // restore
3329    memcpy(solution_,saveSolution,
3330           (numberRows_+numberColumns_)*sizeof(double));
3331    memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char));
3332    memcpy(lower_,saveLower,
3333           (numberRows_+numberColumns_)*sizeof(double));
3334    memcpy(upper_,saveUpper,
3335           (numberRows_+numberColumns_)*sizeof(double));
3336    memcpy(cost_,saveObjective,
3337         (numberRows_+numberColumns_)*sizeof(double));
3338    columnUpper_[iColumn] = saveBound;
3339    memcpy(pivotVariable_, savePivot, numberRows_*sizeof(int));
3340    delete factorization_;
3341    factorization_ = new ClpFactorization(saveFactorization);
3342
3343    if (status||(problemStatus_==0&&!isDualObjectiveLimitReached())) {
3344      objectiveChange = objectiveValue_-saveObjectiveValue;
3345    } else {
3346      objectiveChange = 1.0e100;
3347    }
3348    newUpper[i]=objectiveChange;
3349#ifdef CLP_DEBUG
3350    printf("down on %d costs %g\n",iColumn,objectiveChange);
3351#endif
3352         
3353    // try up
3354   
3355    saveBound = columnLower_[iColumn];
3356    // external view - in case really getting optimal
3357    columnLower_[iColumn] = newLower[i];
3358    if (scalingFlag_<=0) 
3359      lower_[iColumn] = newLower[i];
3360    else 
3361      lower_[iColumn] = newLower[i]/columnScale_[iColumn]; // scale
3362    // Start of fast iterations
3363    status = fastDual(alwaysFinish);
3364    if (status) {
3365      // not finished - might be optimal
3366      checkPrimalSolution(rowActivityWork_,columnActivityWork_);
3367      double limit = 0.0;
3368      getDblParam(ClpDualObjectiveLimit, limit);
3369      if (!numberPrimalInfeasibilities_&&objectiveValue()*optimizationDirection_<
3370          optimizationDirection_*limit) { 
3371        problemStatus_=0;
3372      } 
3373      status=problemStatus_;
3374    }
3375    if (scalingFlag_<=0) {
3376      memcpy(outputSolution[iSolution],solution_,numberColumns_*sizeof(double));
3377    } else {
3378      int j;
3379      double * sol = outputSolution[iSolution];
3380      for (j=0;j<numberColumns_;j++) 
3381        sol[j] = solution_[j]*columnScale_[j];
3382    }
3383    outputStatus[iSolution]=status;
3384    outputIterations[iSolution]=numberIterations_;
3385    iSolution++;
3386
3387    // restore
3388    memcpy(solution_,saveSolution,
3389           (numberRows_+numberColumns_)*sizeof(double));
3390    memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char));
3391    memcpy(lower_,saveLower,
3392           (numberRows_+numberColumns_)*sizeof(double));
3393    memcpy(upper_,saveUpper,
3394           (numberRows_+numberColumns_)*sizeof(double));
3395    memcpy(cost_,saveObjective,
3396         (numberRows_+numberColumns_)*sizeof(double));
3397    columnLower_[iColumn] = saveBound;
3398    memcpy(pivotVariable_, savePivot, numberRows_*sizeof(int));
3399    delete factorization_;
3400    factorization_ = new ClpFactorization(saveFactorization);
3401
3402    if (status||(problemStatus_==0&&!isDualObjectiveLimitReached())) {
3403      objectiveChange = objectiveValue_-saveObjectiveValue;
3404    } else {
3405      objectiveChange = 1.0e100;
3406    }
3407    newLower[i]=objectiveChange;
3408#ifdef CLP_DEBUG
3409    printf("up on %d costs %g\n",iColumn,objectiveChange);
3410#endif
3411         
3412    /* Possibilities are:
3413       Both sides feasible - store
3414       Neither side feasible - set objective high and exit if desired
3415       One side feasible - change bounds and resolve
3416    */
3417    if (newUpper[i]<1.0e100) {
3418      if(newLower[i]<1.0e100) {
3419        // feasible - no action
3420      } else {
3421        // up feasible, down infeasible
3422        returnCode=1;
3423        if (stopOnFirstInfeasible)
3424          break;
3425      }
3426    } else {
3427      if(newLower[i]<1.0e100) {
3428        // down feasible, up infeasible
3429        returnCode=1;
3430        if (stopOnFirstInfeasible)
3431          break;
3432      } else {
3433        // neither side feasible
3434        returnCode=-1;
3435        break;
3436      }
3437    }
3438  }
3439  delete [] saveSolution;
3440  delete [] saveLower;
3441  delete [] saveUpper;
3442  delete [] saveObjective;
3443  delete [] saveStatus;
3444  delete [] savePivot;
3445
3446  // Get rid of some arrays and empty factorization
3447  deleteRim();
3448#else
3449  // save basis and solution
3450  double * rowSolution = new double[numberRows_];
3451  memcpy(rowSolution,rowActivity_,numberRows_*sizeof(double));
3452  double * columnSolution = new double[numberColumns_];
3453  memcpy(columnSolution,columnActivity_,numberColumns_*sizeof(double));
3454  unsigned char * saveStatus = 
3455    new unsigned char [numberRows_+numberColumns_];
3456  if (!status_) {
3457    // odd but anyway setup all slack basis
3458    status_ = new unsigned char [numberColumns_+numberRows_];
3459    memset(status_,0,(numberColumns_+numberRows_)*sizeof(char));
3460  }
3461  memcpy(saveStatus,status_,(numberColumns_+numberRows_)*sizeof(char));
3462  int iSolution =0;
3463  for (i=0;i<numberVariables;i++) {
3464    int iColumn = variables[i];
3465    double objectiveChange;
3466   
3467    // try down
3468   
3469    double saveUpper = columnUpper_[iColumn];
3470    columnUpper_[iColumn] = newUpper[i];
3471    int status=dual(0);
3472    memcpy(outputSolution[iSolution],columnActivity_,numberColumns_*sizeof(double));
3473    outputStatus[iSolution]=status;
3474    outputIterations[iSolution]=numberIterations_;
3475    iSolution++;
3476
3477    // restore
3478    columnUpper_[iColumn] = saveUpper;
3479    memcpy(rowActivity_,rowSolution,numberRows_*sizeof(double));
3480    memcpy(columnActivity_,columnSolution,numberColumns_*sizeof(double));
3481    memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char));
3482
3483    if (problemStatus_==0&&!isDualObjectiveLimitReached()) {
3484      objectiveChange = objectiveValue_-saveObjectiveValue;
3485    } else {
3486      objectiveChange = 1.0e100;
3487    }
3488    newUpper[i]=objectiveChange;
3489#ifdef CLP_DEBUG
3490    printf("down on %d costs %g\n",iColumn,objectiveChange);
3491#endif
3492         
3493    // try up
3494   
3495    double saveLower = columnLower_[iColumn];
3496    columnLower_[iColumn] = newLower[i];
3497    status=dual(0);
3498    memcpy(outputSolution[iSolution],columnActivity_,numberColumns_*sizeof(double));
3499    outputStatus[iSolution]=status;
3500    outputIterations[iSolution]=numberIterations_;
3501    iSolution++;
3502
3503    // restore
3504    columnLower_[iColumn] = saveLower;
3505    memcpy(rowActivity_,rowSolution,numberRows_*sizeof(double));
3506    memcpy(columnActivity_,columnSolution,numberColumns_*sizeof(double));
3507    memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char));
3508
3509    if (problemStatus_==0&&!isDualObjectiveLimitReached()) {
3510      objectiveChange = objectiveValue_-saveObjectiveValue;
3511    } else {
3512      objectiveChange = 1.0e100;
3513    }
3514    newLower[i]=objectiveChange;
3515#ifdef CLP_DEBUG
3516    printf("up on %d costs %g\n",iColumn,objectiveChange);
3517#endif
3518         
3519    /* Possibilities are:
3520       Both sides feasible - store
3521       Neither side feasible - set objective high and exit
3522       One side feasible - change bounds and resolve
3523    */
3524    if (newUpper[i]<1.0e100) {
3525      if(newLower[i]<1.0e100) {
3526        // feasible - no action
3527      } else {
3528        // up feasible, down infeasible
3529        returnCode=1;
3530        if (stopOnFirstInfeasible)
3531          break;
3532      }
3533    } else {
3534      if(newLower[i]<1.0e100) {
3535        // down feasible, up infeasible
3536        returnCode=1;
3537        if (stopOnFirstInfeasible)
3538          break;
3539      } else {
3540        // neither side feasible
3541        returnCode=-1;
3542        break;
3543      }
3544    }
3545  }
3546  delete [] rowSolution;
3547  delete [] columnSolution;
3548  delete [] saveStatus;
3549#endif
3550  objectiveValue_ = saveObjectiveValue;
3551  return returnCode;
3552}
3553// treat no pivot as finished (unless interesting)
3554int ClpSimplexDual::fastDual(bool alwaysFinish)
3555{
3556  algorithm_ = -1;
3557  // save data
3558  ClpDataSave data = saveData();
3559  dualTolerance_=dblParam_[ClpDualTolerance];
3560  primalTolerance_=dblParam_[ClpPrimalTolerance];
3561
3562  // save dual bound
3563  double saveDualBound = dualBound_;
3564
3565  double objectiveChange;
3566  // for dual we will change bounds using dualBound_
3567  // for this we need clean basis so it is after factorize
3568  gutsOfSolution(NULL,NULL);
3569  numberFake_ =0; // Number of variables at fake bounds
3570  changeBounds(true,NULL,objectiveChange);
3571
3572  problemStatus_ = -1;
3573  numberIterations_=0;
3574  factorization_->sparseThreshold(0);
3575  factorization_->goSparse();
3576
3577  int lastCleaned=0; // last time objective or bounds cleaned up
3578
3579  // number of times we have declared optimality
3580  numberTimesOptimal_=0;
3581
3582  // This says whether to restore things etc
3583  int factorType=0;
3584  /*
3585    Status of problem:
3586    0 - optimal
3587    1 - infeasible
3588    2 - unbounded
3589    -1 - iterating
3590    -2 - factorization wanted
3591    -3 - redo checking without factorization
3592    -4 - looks infeasible
3593
3594    BUT also from whileIterating return code is:
3595
3596   -1 iterations etc
3597   -2 inaccuracy
3598   -3 slight inaccuracy (and done iterations)
3599   +0 looks optimal (might be unbounded - but we will investigate)
3600   +1 looks infeasible
3601   +3 max iterations
3602
3603  */
3604
3605  int returnCode = 0;
3606
3607  int iRow,iColumn;
3608  while (problemStatus_<0) {
3609    // clear
3610    for (iRow=0;iRow<4;iRow++) {
3611      rowArray_[iRow]->clear();
3612    }   
3613   
3614    for (iColumn=0;iColumn<2;iColumn++) {
3615      columnArray_[iColumn]->clear();
3616    }   
3617
3618    // give matrix (and model costs and bounds a chance to be
3619    // refreshed (normally null)
3620    matrix_->refresh(this);
3621    // may factorize, checks if problem finished
3622    // should be able to speed this up on first time
3623    statusOfProblemInDual(lastCleaned,factorType,progress_,NULL);
3624
3625    // Say good factorization
3626    factorType=1;
3627
3628    // Do iterations
3629    if (problemStatus_<0) {
3630      double * givenPi=NULL;
3631      returnCode = whileIterating(givenPi);
3632      if (!alwaysFinish&&(returnCode<1||returnCode==3)) {
3633        double limit = 0.0;
3634        getDblParam(ClpDualObjectiveLimit, limit);
3635        if(fabs(limit)>1.0e30||objectiveValue()*optimizationDirection_<
3636           optimizationDirection_*limit|| 
3637           numberAtFakeBound()) {
3638          returnCode=1;
3639          secondaryStatus_ = 1; // and say was on cutoff
3640          // can't say anything interesting - might as well return
3641#ifdef CLP_DEBUG
3642          printf("returning from fastDual after %d iterations with code %d\n",
3643                 numberIterations_,returnCode);
3644#endif
3645          break;
3646        }
3647      }
3648      returnCode=0;
3649    }
3650  }
3651
3652  // clear
3653  for (iRow=0;iRow<4;iRow++) {
3654    rowArray_[iRow]->clear();
3655  }   
3656 
3657  for (iColumn=0;iColumn<2;iColumn++) {
3658    columnArray_[iColumn]->clear();
3659  }   
3660  assert(!numberFake_||returnCode||problemStatus_); // all bounds should be okay
3661  // Restore any saved stuff
3662  restoreData(data);
3663  dualBound_ = saveDualBound;
3664  return returnCode;
3665}
3666/* Checks number of variables at fake bounds.  This is used by fastDual
3667   so can exit gracefully before end */
3668int 
3669ClpSimplexDual::numberAtFakeBound()
3670{
3671  int iSequence;
3672  int numberFake=0;
3673 
3674  for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
3675    FakeBound bound = getFakeBound(iSequence);
3676    switch(getStatus(iSequence)) {
3677
3678    case basic:
3679      break;
3680    case isFree:
3681    case superBasic:
3682    case ClpSimplex::isFixed:
3683      assert (bound==noFake);
3684      break;
3685    case atUpperBound:
3686      if (bound==upperFake||bound==bothFake)
3687        numberFake++;
3688      break;
3689    case atLowerBound:
3690      if (bound==lowerFake||bound==bothFake)
3691        numberFake++;
3692      break;
3693    }
3694  }
3695  numberFake_ = numberFake;
3696  return numberFake;
3697}
3698/* Pivot out a variable and choose an incoing one.  Assumes dual
3699   feasible - will not go through a reduced cost. 
3700   Returns step length in theta
3701   Returns ray in ray_ (or NULL if no pivot)
3702   Return codes as before but -1 means no acceptable pivot
3703*/
3704int 
3705ClpSimplexDual::pivotResult()
3706{
3707  abort();
3708  return 0;
3709}
3710/*
3711   Row array has row part of pivot row
3712   Column array has column part.
3713   This is used in dual values pass
3714*/
3715int
3716ClpSimplexDual::checkPossibleValuesMove(CoinIndexedVector * rowArray,
3717                                        CoinIndexedVector * columnArray,
3718                                        double acceptablePivot,
3719                                        CoinBigIndex * dubiousWeights)
3720{
3721  double * work;
3722  int number;
3723  int * which;
3724  int iSection;
3725
3726  double tolerance = dualTolerance_*1.001;
3727
3728  double thetaDown = 1.0e31;
3729  double changeDown ;
3730  double thetaUp = 1.0e31;
3731  double bestAlphaDown = acceptablePivot*0.99999;
3732  double bestAlphaUp = acceptablePivot*0.99999;
3733  int sequenceDown =-1;
3734  int sequenceUp = sequenceOut_;
3735
3736  double djBasic = dj_[sequenceOut_];
3737  if (djBasic>0.0) {
3738    // basic at lower bound so directionOut_ 1 and -1 in pivot row
3739    // dj will go to zero on other way
3740    thetaUp = djBasic;
3741    changeDown = -lower_[sequenceOut_];
3742  } else {
3743    // basic at upper bound so directionOut_ -1 and 1 in pivot row
3744    // dj will go to zero on other way
3745    thetaUp = -djBasic;
3746    changeDown = upper_[sequenceOut_];
3747  }
3748  bestAlphaUp = 1.0;
3749  int addSequence;
3750
3751  double alphaUp=0.0;
3752  double alphaDown=0.0;
3753
3754  for (iSection=0;iSection<2;iSection++) {
3755
3756    int i;
3757    if (!iSection) {
3758      work = rowArray->denseVector();
3759      number = rowArray->getNumElements();
3760      which = rowArray->getIndices();
3761      addSequence = numberColumns_;
3762    } else {
3763      work = columnArray->denseVector();
3764      number = columnArray->getNumElements();
3765      which = columnArray->getIndices();
3766      addSequence = 0;
3767    }
3768   
3769    for (i=0;i<number;i++) {
3770      int iSequence = which[i];
3771      int iSequence2 = iSequence + addSequence;
3772      double alpha;
3773      double oldValue;
3774      double value;
3775
3776      switch(getStatus(iSequence2)) {
3777         
3778      case basic:
3779        break;
3780      case ClpSimplex::isFixed:
3781        alpha = work[i];
3782        changeDown += alpha*upper_[iSequence2];
3783        break;
3784      case isFree:
3785      case superBasic:
3786        alpha = work[i];
3787        // dj must be effectively zero as dual feasible
3788        if (fabs(alpha)>bestAlphaUp) {
3789          thetaDown = 0.0;
3790          thetaUp = 0.0;
3791          bestAlphaDown = fabs(alpha);
3792          bestAlphaUp = bestAlphaUp;
3793          sequenceDown =iSequence2;
3794          sequenceUp = sequenceDown;
3795          alphaUp = alpha;
3796          alphaDown = alpha;
3797        }
3798        break;
3799      case atUpperBound:
3800        alpha = work[i];
3801        oldValue = dj_[iSequence2];
3802        changeDown += alpha*upper_[iSequence2];
3803        if (alpha>=acceptablePivot) {
3804          // might do other way
3805          value = oldValue+thetaUp*alpha;
3806          if (value>-tolerance) {
3807            if (value>tolerance||fabs(alpha)>bestAlphaUp) {
3808              thetaUp = -oldValue/alpha;
3809              bestAlphaUp = fabs(alpha);
3810              sequenceUp = iSequence2;
3811              alphaUp = alpha;
3812            }
3813          }
3814        } else if (alpha<=-acceptablePivot) {
3815          // might do this way
3816          value = oldValue-thetaDown*alpha;
3817          if (value>-tolerance) {
3818            if (value>tolerance||fabs(alpha)>bestAlphaDown) {
3819              thetaDown = oldValue/alpha;
3820              bestAlphaDown = fabs(alpha);
3821              sequenceDown = iSequence2;
3822              alphaDown = alpha;
3823            }
3824          }
3825        }
3826        break;
3827      case atLowerBound:
3828        alpha = work[i];
3829        oldValue = dj_[iSequence2];
3830        changeDown += alpha*lower_[iSequence2];
3831        if (alpha<=-acceptablePivot) {
3832          // might do other way
3833          value = oldValue+thetaUp*alpha;
3834          if (value<tolerance) {
3835            if (value<-tolerance||fabs(alpha)>bestAlphaUp) {
3836              thetaUp = -oldValue/alpha;
3837              bestAlphaUp = fabs(alpha);
3838              sequenceUp = iSequence2;
3839              alphaUp = alpha;
3840            }
3841          }
3842        } else if (alpha>=acceptablePivot) {
3843          // might do this way
3844          value = oldValue-thetaDown*alpha;
3845          if (value<tolerance) {
3846            if (value<-tolerance||fabs(alpha)>bestAlphaDown) {
3847              thetaDown = oldValue/alpha;
3848              bestAlphaDown = fabs(alpha);
3849              sequenceDown = iSequence2;
3850              alphaDown = alpha;
3851            }
3852          }
3853        }
3854        break;
3855      }
3856    }
3857  }
3858  int returnCode;
3859  thetaUp *= -1.0;
3860  double changeUp = -thetaUp * changeDown;
3861  changeDown = -thetaDown*changeDown;
3862  if (max(fabs(thetaDown),fabs(thetaUp))<1.0e-8) {
3863    // largest
3864    if (fabs(alphaDown)<fabs(alphaUp)) {
3865      sequenceDown =-1;
3866    }
3867  }
3868  // choose
3869  if (changeDown>changeUp&&sequenceDown>=0) {
3870    theta_ = thetaDown;
3871    sequenceIn_ = sequenceDown;
3872    alpha_ = alphaDown;
3873#ifdef CLP_DEBUG
3874    if ((handler_->logLevel()&32))
3875      printf("predicted way - dirout %d, change %g,%g theta %g\n",
3876             directionOut_,changeDown,changeUp,theta_);
3877#endif
3878    returnCode = 0;
3879  } else {
3880    theta_ = thetaUp;
3881    sequenceIn_ = sequenceUp;
3882    alpha_ = alphaUp;
3883    if (sequenceIn_!=sequenceOut_) {
3884#ifdef CLP_DEBUG
3885      if ((handler_->logLevel()&32))
3886        printf("opposite way - dirout %d, change %g,%g theta %g\n",
3887               directionOut_,changeDown,changeUp,theta_);
3888#endif
3889      returnCode = 0;
3890    } else {
3891#ifdef CLP_DEBUG
3892      if ((handler_->logLevel()&32))
3893        printf("opposite way to zero dj - dirout %d, change %g,%g theta %g\n",
3894               directionOut_,changeDown,changeUp,theta_);
3895#endif
3896      returnCode = 1;
3897    }
3898  }
3899  return returnCode;
3900}
3901/*
3902   This sees if we can move duals in dual values pass.
3903   This is done before any pivoting
3904*/
3905void ClpSimplexDual::doEasyOnesInValuesPass(double * dj)
3906{
3907  // Get column copy
3908  CoinPackedMatrix * columnCopy = matrix();
3909  // Get a row copy in standard format
3910  CoinPackedMatrix copy;
3911  copy.reverseOrderedCopyOf(*columnCopy);
3912  // get matrix data pointers
3913  const int * column = copy.getIndices();
3914  const CoinBigIndex * rowStart = copy.getVectorStarts();
3915  const int * rowLength = copy.getVectorLengths(); 
3916  const double * elementByRow = copy.getElements();
3917  double tolerance = dualTolerance_*1.001;
3918
3919  int iRow;
3920#ifdef CLP_DEBUG
3921  {
3922    double value5=0.0;
3923    int i;
3924    for (i=0;i<numberRows_+numberColumns_;i++) {
3925      if (dj[i]<-1.0e-6)
3926        value5 += dj[i]*upper_[i];
3927      else if (dj[i] >1.0e-6)
3928        value5 += dj[i]*lower_[i];
3929    }
3930    printf("Values objective Value before %g\n",value5);
3931  }
3932#endif
3933  // for scaled row
3934  double * scaled=NULL;
3935  if (rowScale_)
3936    scaled = new double[numberColumns_];
3937  for (iRow=0;iRow<numberRows_;iRow++) {
3938
3939    int iSequence = iRow + numberColumns_;
3940    double djBasic = dj[iSequence];
3941    if (getRowStatus(iRow)==basic&&fabs(djBasic)>tolerance) {
3942
3943      double changeUp ;
3944      // always -1 in pivot row
3945      if (djBasic>0.0) {
3946        // basic at lower bound
3947        changeUp = -lower_[iSequence];
3948      } else {
3949        // basic at upper bound
3950        changeUp = upper_[iSequence];
3951      }
3952      bool canMove=true;
3953      int i;
3954      const double * thisElements = elementByRow + rowStart[iRow]; 
3955      const int * thisIndices = column+rowStart[iRow];
3956      if (rowScale_) {
3957        // scale row
3958        double scale = rowScale_[iRow];
3959        for (i = 0;i<rowLength[iRow];i++) {
3960          int iColumn = thisIndices[i];
3961          double alpha = thisElements[i];
3962          scaled[i] = scale*alpha*columnScale_[iColumn];
3963        }
3964        thisElements = scaled;
3965      }
3966      for (i = 0;i<rowLength[iRow];i++) {
3967        int iColumn = thisIndices[i];
3968        double alpha = thisElements[i];
3969        double oldValue = dj[iColumn];;
3970        double value;
3971
3972        switch(getStatus(iColumn)) {
3973         
3974        case basic:
3975          if (dj[iColumn]<-tolerance&&
3976              fabs(solution_[iColumn]-upper_[iColumn])<1.0e-8) {
3977            // at ub
3978            changeUp += alpha*upper_[iColumn];
3979            // might do other way
3980            value = oldValue+djBasic*alpha;
3981            if (value>tolerance) 
3982              canMove=false;
3983          } else if (dj[iColumn]>tolerance&&
3984              fabs(solution_[iColumn]-lower_[iColumn])<1.0e-8) {
3985            changeUp += alpha*lower_[iColumn];
3986            // might do other way
3987            value = oldValue+djBasic*alpha;
3988            if (value<-tolerance)
3989              canMove=false;
3990          } else {
3991            canMove=false;
3992          }
3993          break;
3994        case ClpSimplex::isFixed:
3995          changeUp += alpha*upper_[iColumn];
3996          break;
3997        case isFree:
3998        case superBasic:
3999          canMove=false;
4000        break;
4001      case atUpperBound:
4002        changeUp += alpha*upper_[iColumn];
4003        // might do other way
4004        value = oldValue+djBasic*alpha;
4005        if (value>tolerance) 
4006          canMove=false;
4007        break;
4008      case atLowerBound:
4009        changeUp += alpha*lower_[iColumn];
4010        // might do other way
4011        value = oldValue+djBasic*alpha;
4012        if (value<-tolerance)
4013          canMove=false;
4014        break;
4015        }
4016      }
4017      if (canMove) {
4018        if (changeUp*djBasic>1.0e-12||fabs(changeUp)<1.0e-8) {
4019          // move
4020          for (i = 0;i<rowLength[iRow];i++) {
4021            int iColumn = thisIndices[i];
4022            double alpha = thisElements[i];
4023            dj[iColumn] += djBasic * alpha;
4024          }
4025          dj[iSequence]=0.0;
4026#ifdef CLP_DEBUG
4027          {
4028            double value5=0.0;
4029            int i;
4030            for (i=0;i<numberRows_+numberColumns_;i++) {
4031              if (dj[i]<-1.0e-6)
4032                value5 += dj[i]*upper_[i];
4033              else if (dj[i] >1.0e-6)
4034                value5 += dj[i]*lower_[i];
4035            }
4036            printf("Values objective Value after row %d old dj %g %g\n",
4037                   iRow,djBasic,value5);
4038          }
4039#endif
4040        }
4041      }
4042    }     
4043  }
4044  delete [] scaled;
4045}
4046int
4047ClpSimplexDual::nextSuperBasic()
4048{
4049  if (firstFree_>=0) {
4050    int returnValue=firstFree_;
4051    int iColumn=firstFree_+1;
4052    for (;iColumn<numberRows_+numberColumns_;iColumn++) {
4053      if (getStatus(iColumn)==isFree) 
4054        if (fabs(dj_[iColumn])>1.0e2*dualTolerance_) 
4055          break;
4056    }
4057    firstFree_ = iColumn;
4058    if (firstFree_==numberRows_+numberColumns_)
4059      firstFree_=-1;
4060    return returnValue;
4061  } else {
4062    return -1;
4063  }
4064}
Note: See TracBrowser for help on using the repository browser.