source: trunk/ClpSimplexDual.cpp @ 403

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

Bug - Francois margot

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