source: branches/pre/ClpSimplexDual.cpp @ 210

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

Trying

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