source: trunk/ClpSimplexDual.cpp @ 361

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

Saving numberFake_

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