source: trunk/ClpSimplexDual.cpp @ 348

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

New status 5

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