source: trunk/ClpSimplexDual.cpp @ 437

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

secondary status

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