source: trunk/ClpSimplexDual.cpp @ 553

Last change on this file since 553 was 553, checked in by forrest, 15 years ago

try and make dual infeas faster

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