# source:trunk/Clp/src/ClpSimplexDual.cpp@1277

Last change on this file since 1277 was 1277, checked in by forrest, 11 years ago

• Property svn:eol-style set to `native`
• Property svn:keywords set to `Author Date Id Revision`
File size: 189.1 KB
Line
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
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#ifdef COIN_DEVELOP
89#undef COIN_DEVELOP
90#define COIN_DEVELOP 2
91#endif
92
93#include "CoinPragma.hpp"
94
95#include <math.h>
96
97#include "CoinHelperFunctions.hpp"
98#include "ClpHelperFunctions.hpp"
99#include "ClpSimplexDual.hpp"
100#include "ClpEventHandler.hpp"
101#include "ClpFactorization.hpp"
102#include "CoinPackedMatrix.hpp"
103#include "CoinIndexedVector.hpp"
104#include "ClpDualRowDantzig.hpp"
105#include "ClpMessage.hpp"
106#include "ClpLinearObjective.hpp"
107#include <cfloat>
108#include <cassert>
109#include <string>
110#include <stdio.h>
111#include <iostream>
112//#define CLP_DEBUG 1
113// To force to follow another run put logfile name here and define
114//#define FORCE_FOLLOW
115#ifdef FORCE_FOLLOW
116static FILE * fpFollow=NULL;
117static char * forceFile="old.log";
118static int force_in=-1;
119static int force_out=-1;
120static int force_iteration=0;
121#endif
122//#define VUB
123#ifdef VUB
124extern int * vub;
125extern int * toVub;
126extern int * nextDescendent;
127#endif
128// dual
129
130  /* *** Method
131     This is a vanilla version of dual simplex.
132
133     It tries to be a single phase approach with a weight of 1.0 being
134     given to getting optimal and a weight of dualBound_ being
135     given to getting dual feasible.  In this version I have used the
136     idea that this weight can be thought of as a fake bound.  If the
137     distance between the lower and upper bounds on a variable is less
138     than the feasibility weight then we are always better off flipping
139     to other bound to make dual feasible.  If the distance is greater
140     then we make up a fake bound dualBound_ away from one bound.
141     If we end up optimal or primal infeasible, we check to see if
142     bounds okay.  If so we have finished, if not we increase dualBound_
143     and continue (after checking if unbounded). I am undecided about
144     free variables - there is coding but I am not sure about it.  At
145     present I put them in basis anyway.
146
147     The code is designed to take advantage of sparsity so arrays are
148     seldom zeroed out from scratch or gone over in their entirety.
149     The only exception is a full scan to find outgoing variable.  This
150     will be changed to keep an updated list of infeasibilities (or squares
151     if steepest edge).  Also on easy problems we don't need full scan - just
152     pick first reasonable.
153
154     One problem is how to tackle degeneracy and accuracy.  At present
155     I am using the modification of costs which I put in OSL and which was
156     extended by Gill et al.  I am still not sure of the exact details.
157
158     The flow of dual is three while loops as follows:
159
160     while (not finished) {
161
162       while (not clean solution) {
163
164          Factorize and/or clean up solution by flipping variables so
165          dual feasible.  If looks finished check fake dual bounds.
166          Repeat until status is iterating (-1) or finished (0,1,2)
167
168       }
169
170       while (status==-1) {
171
172         Iterate until no pivot in or out or time to re-factorize.
173
174         Flow is:
175
176         choose pivot row (outgoing variable).  if none then
177         we are primal feasible so looks as if done but we need to
178         break and check bounds etc.
179
180         Get pivot row in tableau
181
182         Choose incoming column.  If we don't find one then we look
183         primal infeasible so break and check bounds etc.  (Also the
184         pivot tolerance is larger after any iterations so that may be
185         reason)
186
187         If we do find incoming column, we may have to adjust costs to
188         keep going forwards (anti-degeneracy).  Check pivot will be stable
189         and if unstable throw away iteration (we will need to implement
190         flagging of basic variables sometime) and break to re-factorize.
191         If minor error re-factorize after iteration.
192
193         Update everything (this may involve flipping variables to stay
194         dual feasible.
195
196       }
197
198     }
199
200     At present we never check we are going forwards.  I overdid that in
201     OSL so will try and make a last resort.
202
203     Needs partial scan pivot out option.
204     Needs dantzig, uninitialized and full steepest edge options (can still
205     use partial scan)
206
207     May need other anti-degeneracy measures, especially if we try and use
208     loose tolerances as a way to solve in fewer iterations.
209
210     I like idea of dynamic scaling.  This gives opportunity to decouple
211     different implications of scaling for accuracy, iteration count and
212     feasibility tolerance.
213
214  */
215#define CLEAN_FIXED 0
216// Startup part of dual (may be extended to other algorithms)
217int
218ClpSimplexDual::startupSolve(int ifValuesPass,double * saveDuals,int startFinishOptions)
219{
220  // If values pass then save given duals round check solution
221  // sanity check
222  // initialize - no values pass and algorithm_ is -1
223  // put in standard form (and make row copy)
224  // create modifiable copies of model rim and do optional scaling
225  // If problem looks okay
226  // Do initial factorization
227  // If user asked for perturbation - do it
228  if (!startup(0,startFinishOptions)) {
229    int usePrimal=0;
230    // looks okay
231    // Superbasic variables not allowed
232    // If values pass then scale pi
233    if (ifValuesPass) {
234      if (problemStatus_&&perturbation_<100)
235        usePrimal=perturb();
236      int i;
237      if (scalingFlag_>0) {
238        for (i=0;i<numberRows_;i++) {
239          dual_[i] = saveDuals[i]/rowScale_[i];
240        }
241      } else {
242        CoinMemcpyN(saveDuals,numberRows_,dual_);
243      }
244      // now create my duals
245      for (i=0;i<numberRows_;i++) {
246        // slack
247        double value = dual_[i];
248        value += rowObjectiveWork_[i];
249        saveDuals[i+numberColumns_]=value;
250      }
251      CoinMemcpyN(objectiveWork_,numberColumns_,saveDuals);
252      transposeTimes(-1.0,dual_,saveDuals);
253      // make reduced costs okay
254      for (i=0;i<numberColumns_;i++) {
255        if (getStatus(i)==atLowerBound) {
256          if (saveDuals[i]<0.0) {
257            //if (saveDuals[i]<-1.0e-3)
258            //printf("bad dj at lb %d %g\n",i,saveDuals[i]);
259            saveDuals[i]=0.0;
260          }
261        } else if (getStatus(i)==atUpperBound) {
262          if (saveDuals[i]>0.0) {
263            //if (saveDuals[i]>1.0e-3)
264            //printf("bad dj at ub %d %g\n",i,saveDuals[i]);
265            saveDuals[i]=0.0;
266          }
267        }
268      }
269      CoinMemcpyN(saveDuals,(numberColumns_+numberRows_),dj_);
270      // set up possible ones
271      for (i=0;i<numberRows_+numberColumns_;i++)
272        clearPivoted(i);
273      int iRow;
274      for (iRow=0;iRow<numberRows_;iRow++) {
275        int iPivot=pivotVariable_[iRow];
276        if (fabs(saveDuals[iPivot])>dualTolerance_) {
277          if (getStatus(iPivot)!=isFree)
278            setPivoted(iPivot);
279        }
280      }
281    } else if ((specialOptions_&1024)!=0&&CLEAN_FIXED) {
282      // set up possible ones
283      for (int i=0;i<numberRows_+numberColumns_;i++)
284        clearPivoted(i);
285      int iRow;
286      for (iRow=0;iRow<numberRows_;iRow++) {
287        int iPivot=pivotVariable_[iRow];
288        if (iPivot<numberColumns_&&lower_[iPivot]==upper_[iPivot]) {
289          setPivoted(iPivot);
290        }
291      }
292    }
293
294    double objectiveChange;
295    numberFake_ =0; // Number of variables at fake bounds
296    numberChanged_ =0; // Number of variables with changed costs
297    changeBounds(1,NULL,objectiveChange);
298
299    if (!ifValuesPass) {
300      // Check optimal
301      if (!numberDualInfeasibilities_&&!numberPrimalInfeasibilities_)
302        problemStatus_=0;
303    }
304    if (problemStatus_<0&&perturbation_<100) {
305      bool inCbcOrOther = (specialOptions_&0x03000000)!=0;
306      if (!inCbcOrOther)
307        usePrimal=perturb();
308      // Can't get here if values pass
309      gutsOfSolution(NULL,NULL);
310      if (handler_->logLevel()>2) {
311        handler_->message(CLP_SIMPLEX_STATUS,messages_)
312          <<numberIterations_<<objectiveValue();
313        handler_->printing(sumPrimalInfeasibilities_>0.0)
314          <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
315        handler_->printing(sumDualInfeasibilities_>0.0)
316          <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
317        handler_->printing(numberDualInfeasibilitiesWithoutFree_
318                           <numberDualInfeasibilities_)
319                             <<numberDualInfeasibilitiesWithoutFree_;
320        handler_->message()<<CoinMessageEol;
321      }
322      if (inCbcOrOther) {
323        if (numberPrimalInfeasibilities_) {
324          usePrimal=perturb();
325          if (perturbation_>=101)
326            computeDuals(NULL);
327          //gutsOfSolution(NULL,NULL);
328        } else if (numberDualInfeasibilities_) {
329          problemStatus_=10;
330#if COIN_DEVELOP>2
331
332          printf("returning at %d\n",__LINE__);
333#endif
334          return 1; // to primal
335        }
336      }
337    } else if (!ifValuesPass) {
338      gutsOfSolution(NULL,NULL);
339    }
340    if (usePrimal) {
341      problemStatus_=10;
342#if COIN_DEVELOP>1
343      printf("returning to use primal (no obj) at %d\n",__LINE__);
344#endif
345    }
346    return usePrimal;
347  } else {
348    return 1;
349  }
350}
351void
352ClpSimplexDual::finishSolve(int startFinishOptions)
353{
354  assert (problemStatus_||!sumPrimalInfeasibilities_);
355
356  // clean up
357  finish(startFinishOptions);
358}
359#ifdef CLP_INVESTIGATE
360static int z_reason[7]={0,0,0,0,0,0,0};
361static int z_thinks=-1;
362#endif
363void
364ClpSimplexDual::gutsOfDual(int ifValuesPass,double * & saveDuals,int initialStatus,
365                           ClpDataSave & data)
366{
367#ifdef CLP_INVESTIGATE
368  z_reason[0]++;
369  z_thinks=-1;
370  int nPivots=9999;
371#endif
372  int lastCleaned=0; // last time objective or bounds cleaned up
373
374  // This says whether to restore things etc
375  // startup will have factorized so can skip
376  int factorType=0;
377  // Start check for cycles
378  progress_.startCheck();
379  // Say change made on first iteration
381  /*
382    Status of problem:
383    0 - optimal
384    1 - infeasible
385    2 - unbounded
386    -1 - iterating
387    -2 - factorization wanted
388    -3 - redo checking without factorization
389    -4 - looks infeasible
390  */
391  while (problemStatus_<0) {
392    int iRow, iColumn;
393    // clear
394    for (iRow=0;iRow<4;iRow++) {
395      rowArray_[iRow]->clear();
396    }
397
398    for (iColumn=0;iColumn<2;iColumn++) {
399      columnArray_[iColumn]->clear();
400    }
401
402    // give matrix (and model costs and bounds a chance to be
403    // refreshed (normally null)
404    matrix_->refresh(this);
405    // If getting nowhere - why not give it a kick
406    // does not seem to work too well - do some more work
407    if (perturbation_<101&&numberIterations_>2*(numberRows_+numberColumns_)
408        &&initialStatus!=10) {
409      perturb();
410      // Can't get here if values pass
411      gutsOfSolution(NULL,NULL);
412      if (handler_->logLevel()>2) {
413        handler_->message(CLP_SIMPLEX_STATUS,messages_)
414          <<numberIterations_<<objectiveValue();
415        handler_->printing(sumPrimalInfeasibilities_>0.0)
416          <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
417        handler_->printing(sumDualInfeasibilities_>0.0)
418          <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
419        handler_->printing(numberDualInfeasibilitiesWithoutFree_
420                           <numberDualInfeasibilities_)
421                             <<numberDualInfeasibilitiesWithoutFree_;
422        handler_->message()<<CoinMessageEol;
423      }
424    }
425    // see if in Cbc etc
426    bool inCbcOrOther = (specialOptions_&0x03000000)!=0;
427#if 0
428    bool gotoPrimal=false;
429    if (inCbcOrOther&&numberIterations_>disasterArea_+numberRows_&&
430        numberDualInfeasibilitiesWithoutFree_&&largestDualError_>1.0e-1) {
431      if (!disasterArea_) {
432        printf("trying all slack\n");
433        // try all slack basis
434        allSlackBasis(true);
435        disasterArea_=2*numberRows_;
436      } else {
437        printf("going to primal\n");
438        // go to primal
439        gotoPrimal=true;
440        allSlackBasis(true);
441      }
442    }
443#endif
444    bool disaster=false;
445    if (disasterArea_&&inCbcOrOther&&disasterArea_->check()) {
446      disasterArea_->saveInfo();
447      disaster=true;
448    }
449    // may factorize, checks if problem finished
450    statusOfProblemInDual(lastCleaned,factorType,saveDuals,data,
451                          ifValuesPass);
452    if (disaster)
453      problemStatus_=3;
454    // If values pass then do easy ones on first time
455    if (ifValuesPass&&
456        progress_.lastIterationNumber(0)<0&&saveDuals) {
457      doEasyOnesInValuesPass(saveDuals);
458    }
459
460    // Say good factorization
461    factorType=1;
462    if (data.sparseThreshold_) {
463      // use default at present
464      factorization_->sparseThreshold(0);
465      factorization_->goSparse();
466    }
467
468    // exit if victory declared
469    if (problemStatus_>=0)
470      break;
471
472    // test for maximum iterations
473    if (hitMaximumIterations()||(ifValuesPass==2&&!saveDuals)) {
474      problemStatus_=3;
475      break;
476    }
477    if (ifValuesPass&&!saveDuals) {
478      // end of values pass
479      ifValuesPass=0;
480      int status = eventHandler_->event(ClpEventHandler::endOfValuesPass);
481      if (status>=0) {
482        problemStatus_=5;
483        secondaryStatus_=ClpEventHandler::endOfValuesPass;
484        break;
485      }
486    }
487    // Check event
488    {
489      int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
490      if (status>=0) {
491        problemStatus_=5;
492        secondaryStatus_=ClpEventHandler::endOfFactorization;
493        break;
494      }
495    }
496    // Do iterations
497    int returnCode=whileIterating(saveDuals,ifValuesPass);
498#ifdef CLP_INVESTIGATE
499    nPivots=factorization_->pivots();
500#endif
501    if (returnCode==-2)
502      factorType=3;
503  }
504#ifdef CLP_INVESTIGATE
505  if (z_thinks!=-1) {
506    assert (z_thinks<4);
507    if ((!factorization_->pivots()&&nPivots<20)&&z_thinks>=0&&z_thinks<2)
508      z_thinks += 4;
509    z_reason[1+z_thinks]++;
510  }
511  if ((z_reason[0]%1000)==0) {
512    printf("Reason");
513    for (int i=0;i<7;i++)
514      printf(" %d",z_reason[i]);
515    printf("\n");
516  }
517#endif
518}
519int
520ClpSimplexDual::dual(int ifValuesPass,int startFinishOptions)
521{
522  algorithm_ = -1;
523  moreSpecialOptions_ &= ~16; // clear check replaceColumn accuracy
524  // save data
525  ClpDataSave data = saveData();
526  double * saveDuals = NULL;
527  int saveDont = dontFactorizePivots_;
528  if ((specialOptions_&2048)==0)
529    dontFactorizePivots_=0;
530  else if(!dontFactorizePivots_)
531    dontFactorizePivots_ = 20;
532  if (ifValuesPass) {
533    saveDuals = new double [numberRows_+numberColumns_];
534    CoinMemcpyN(dual_,numberRows_,saveDuals);
535  }
536  if (alphaAccuracy_!=-1.0)
537    alphaAccuracy_ = 1.0;
538  int returnCode = startupSolve(ifValuesPass,saveDuals,startFinishOptions);
539  // Save so can see if doing after primal
540  int initialStatus=problemStatus_;
541
542  if (!returnCode)
543    gutsOfDual(ifValuesPass,saveDuals,initialStatus,data);
544  if (problemStatus_==10)
545    startFinishOptions |= 1;
546  finishSolve(startFinishOptions);
547  delete [] saveDuals;
548
549  // Restore any saved stuff
550  restoreData(data);
551  dontFactorizePivots_ = saveDont;
552  return problemStatus_;
553}
554// old way
555#if 0
556int ClpSimplexDual::dual (int ifValuesPass , int startFinishOptions)
557{
558  algorithm_ = -1;
559
560  // save data
561  ClpDataSave data = saveData();
562  // Save so can see if doing after primal
563  int initialStatus=problemStatus_;
564
565  // If values pass then save given duals round check solution
566  double * saveDuals = NULL;
567  if (ifValuesPass) {
568    saveDuals = new double [numberRows_+numberColumns_];
569    CoinMemcpyN(dual_,numberRows_,saveDuals);
570  }
571  // sanity check
572  // initialize - no values pass and algorithm_ is -1
573  // put in standard form (and make row copy)
574  // create modifiable copies of model rim and do optional scaling
575  // If problem looks okay
576  // Do initial factorization
577  // If user asked for perturbation - do it
578  if (!startup(0,startFinishOptions)) {
579    // looks okay
580    // Superbasic variables not allowed
581    // If values pass then scale pi
582    if (ifValuesPass) {
583      if (problemStatus_&&perturbation_<100)
584        perturb();
585      int i;
586      if (scalingFlag_>0) {
587        for (i=0;i<numberRows_;i++) {
588          dual_[i] = saveDuals[i]/rowScale_[i];
589        }
590      } else {
591        CoinMemcpyN(saveDuals,numberRows_,dual_);
592      }
593      // now create my duals
594      for (i=0;i<numberRows_;i++) {
595        // slack
596        double value = dual_[i];
597        value += rowObjectiveWork_[i];
598        saveDuals[i+numberColumns_]=value;
599      }
600      CoinMemcpyN(objectiveWork_,numberColumns_,saveDuals);
601      transposeTimes(-1.0,dual_,saveDuals);
602      // make reduced costs okay
603      for (i=0;i<numberColumns_;i++) {
604        if (getStatus(i)==atLowerBound) {
605          if (saveDuals[i]<0.0) {
606            //if (saveDuals[i]<-1.0e-3)
607            //printf("bad dj at lb %d %g\n",i,saveDuals[i]);
608            saveDuals[i]=0.0;
609          }
610        } else if (getStatus(i)==atUpperBound) {
611          if (saveDuals[i]>0.0) {
612            //if (saveDuals[i]>1.0e-3)
613            //printf("bad dj at ub %d %g\n",i,saveDuals[i]);
614            saveDuals[i]=0.0;
615          }
616        }
617      }
618      CoinMemcpyN(saveDuals,numberColumns_+numberRows_,dj_);
619      // set up possible ones
620      for (i=0;i<numberRows_+numberColumns_;i++)
621        clearPivoted(i);
622      int iRow;
623      for (iRow=0;iRow<numberRows_;iRow++) {
624        int iPivot=pivotVariable_[iRow];
625        if (fabs(saveDuals[iPivot])>dualTolerance_) {
626          if (getStatus(iPivot)!=isFree)
627            setPivoted(iPivot);
628        }
629      }
630    } else if ((specialOptions_&1024)!=0&&CLEAN_FIXED) {
631      // set up possible ones
632      for (int i=0;i<numberRows_+numberColumns_;i++)
633        clearPivoted(i);
634      int iRow;
635      for (iRow=0;iRow<numberRows_;iRow++) {
636        int iPivot=pivotVariable_[iRow];
637        if (iPivot<numberColumns_&&lower_[iPivot]==upper_[iPivot]) {
638          setPivoted(iPivot);
639        }
640      }
641    }
642
643    double objectiveChange;
644    numberFake_ =0; // Number of variables at fake bounds
645    numberChanged_ =0; // Number of variables with changed costs
646    changeBounds(1,NULL,objectiveChange);
647
648    int lastCleaned=0; // last time objective or bounds cleaned up
649
650    if (!ifValuesPass) {
651      // Check optimal
652      if (!numberDualInfeasibilities_&&!numberPrimalInfeasibilities_)
653        problemStatus_=0;
654    }
655    if (problemStatus_<0&&perturbation_<100) {
656      perturb();
657      // Can't get here if values pass
658      gutsOfSolution(NULL,NULL);
659      if (handler_->logLevel()>2) {
660        handler_->message(CLP_SIMPLEX_STATUS,messages_)
661          <<numberIterations_<<objectiveValue();
662        handler_->printing(sumPrimalInfeasibilities_>0.0)
663          <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
664        handler_->printing(sumDualInfeasibilities_>0.0)
665          <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
666        handler_->printing(numberDualInfeasibilitiesWithoutFree_
667                           <numberDualInfeasibilities_)
668                             <<numberDualInfeasibilitiesWithoutFree_;
669        handler_->message()<<CoinMessageEol;
670      }
671    }
672
673    // This says whether to restore things etc
674    // startup will have factorized so can skip
675    int factorType=0;
676    // Start check for cycles
677    progress_.startCheck();
678    // Say change made on first iteration
680    /*
681      Status of problem:
682      0 - optimal
683      1 - infeasible
684      2 - unbounded
685      -1 - iterating
686      -2 - factorization wanted
687      -3 - redo checking without factorization
688      -4 - looks infeasible
689    */
690    while (problemStatus_<0) {
691      int iRow, iColumn;
692      // clear
693      for (iRow=0;iRow<4;iRow++) {
694        rowArray_[iRow]->clear();
695      }
696
697      for (iColumn=0;iColumn<2;iColumn++) {
698        columnArray_[iColumn]->clear();
699      }
700
701      // give matrix (and model costs and bounds a chance to be
702      // refreshed (normally null)
703      matrix_->refresh(this);
704      // If getting nowhere - why not give it a kick
705      // does not seem to work too well - do some more work
706      if (perturbation_<101&&numberIterations_>2*(numberRows_+numberColumns_)
707          &&initialStatus!=10) {
708        perturb();
709        // Can't get here if values pass
710        gutsOfSolution(NULL,NULL);
711        if (handler_->logLevel()>2) {
712          handler_->message(CLP_SIMPLEX_STATUS,messages_)
713            <<numberIterations_<<objectiveValue();
714          handler_->printing(sumPrimalInfeasibilities_>0.0)
715            <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
716          handler_->printing(sumDualInfeasibilities_>0.0)
717            <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
718          handler_->printing(numberDualInfeasibilitiesWithoutFree_
719                             <numberDualInfeasibilities_)
720                               <<numberDualInfeasibilitiesWithoutFree_;
721          handler_->message()<<CoinMessageEol;
722        }
723      }
724      // may factorize, checks if problem finished
725      statusOfProblemInDual(lastCleaned,factorType,saveDuals,data,
726                            ifValuesPass);
727      // If values pass then do easy ones on first time
728      if (ifValuesPass&&
729          progress_.lastIterationNumber(0)<0&&saveDuals) {
730        doEasyOnesInValuesPass(saveDuals);
731      }
732
733      // Say good factorization
734      factorType=1;
735      if (data.sparseThreshold_) {
736        // use default at present
737        factorization_->sparseThreshold(0);
738        factorization_->goSparse();
739      }
740
741      // exit if victory declared
742      if (problemStatus_>=0)
743        break;
744
745      // test for maximum iterations
746      if (hitMaximumIterations()||(ifValuesPass==2&&!saveDuals)) {
747        problemStatus_=3;
748        break;
749      }
750      if (ifValuesPass&&!saveDuals) {
751        // end of values pass
752        ifValuesPass=0;
753        int status = eventHandler_->event(ClpEventHandler::endOfValuesPass);
754        if (status>=0) {
755          problemStatus_=5;
756          secondaryStatus_=ClpEventHandler::endOfValuesPass;
757          break;
758        }
759      }
760      // Check event
761      {
762        int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
763        if (status>=0) {
764          problemStatus_=5;
765          secondaryStatus_=ClpEventHandler::endOfFactorization;
766          break;
767        }
768      }
769      // Do iterations
770      whileIterating(saveDuals,ifValuesPass);
771    }
772  }
773
774  assert (problemStatus_||!sumPrimalInfeasibilities_);
775
776  // clean up
777  finish(startFinishOptions);
778  delete [] saveDuals;
779
780  // Restore any saved stuff
781  restoreData(data);
782  return problemStatus_;
783}
784#endif
785//#define CHECK_ACCURACY
786#ifdef CHECK_ACCURACY
787static double zzzzzz[100000];
788#endif
789/* Reasons to come out:
790   -1 iterations etc
791   -2 inaccuracy
792   -3 slight inaccuracy (and done iterations)
793   +0 looks optimal (might be unbounded - but we will investigate)
794   +1 looks infeasible
795   +3 max iterations
796 */
797int
798ClpSimplexDual::whileIterating(double * & givenDuals,int ifValuesPass)
799{
800#ifdef CLP_INVESTIGATE
801  z_thinks=-1;
802#endif
803#if 0
804  if (!numberIterations_&&auxiliaryModel_) {
805    for (int i=0;i<numberColumns_;i++) {
806      if (!columnLower_[i]==auxiliaryModel_->lowerRegion()[i+numberRows_+numberColumns_])
807        {printf("%d a\n",i);break;}
808      if (!columnUpper_[i]==auxiliaryModel_->upperRegion()[i+numberRows_+numberColumns_])
809        {printf("%d b %g\n",i,auxiliaryModel_->upperRegion()[i+numberRows_+numberColumns_]);break;}
810    }
811    for (int i=0 ;i<numberRows_;i++) {
812      if (!rowLower_[i]==auxiliaryModel_->lowerRegion()[i+numberRows_+2*numberColumns_])
813        {printf("%d c\n",i);break;}
814      if (!rowUpper_[i]==auxiliaryModel_->upperRegion()[i+numberRows_+2*numberColumns_])
815        {printf("%d d\n",i);break;}
816    }
817  }
818#endif
819#ifdef CLP_DEBUG
820  int debugIteration=-1;
821#endif
822  {
823    int i;
824    for (i=0;i<4;i++) {
825      rowArray_[i]->clear();
826    }
827    for (i=0;i<2;i++) {
828      columnArray_[i]->clear();
829    }
830  }
831  // if can't trust much and long way from optimal then relax
832  if (largestPrimalError_>10.0)
833    factorization_->relaxAccuracyCheck(CoinMin(1.0e2,largestPrimalError_/10.0));
834  else
835    factorization_->relaxAccuracyCheck(1.0);
836  // status stays at -1 while iterating, >=0 finished, -2 to invert
837  // status -3 to go to top without an invert
838  int returnCode = -1;
839  double saveSumDual = sumDualInfeasibilities_; // so we know to be careful
840
841#if 0
842  // compute average infeasibility for backward test
843  double averagePrimalInfeasibility = sumPrimalInfeasibilities_/
844    ((double ) (numberPrimalInfeasibilities_+1));
845#endif
846
847  // Get dubious weights
848  CoinBigIndex * dubiousWeights = NULL;
849#ifdef DUBIOUS_WEIGHTS
850  factorization_->getWeights(rowArray_[0]->getIndices());
851  dubiousWeights = matrix_->dubiousWeights(this,rowArray_[0]->getIndices());
852#endif
853  // If values pass then get list of candidates
854  int * candidateList = NULL;
855  int numberCandidates = 0;
856#ifdef CLP_DEBUG
857  bool wasInValuesPass= (givenDuals!=NULL);
858#endif
859  int candidate=-1;
860  if (givenDuals) {
861    assert (ifValuesPass);
862    ifValuesPass=1;
863    candidateList = new int[numberRows_];
864    // move reduced costs across
865    CoinMemcpyN(givenDuals,numberRows_+numberColumns_,dj_);
866    int iRow;
867    for (iRow=0;iRow<numberRows_;iRow++) {
868      int iPivot=pivotVariable_[iRow];
869      if (flagged(iPivot))
870        continue;
871      if (fabs(dj_[iPivot])>dualTolerance_) {
872        // for now safer to ignore free ones
873        if (lower_[iPivot]>-1.0e50||upper_[iPivot]<1.0e50)
874          if (pivoted(iPivot))
875            candidateList[numberCandidates++]=iRow;
876      } else {
877        clearPivoted(iPivot);
878      }
879    }
880    // and set first candidate
881    if (!numberCandidates) {
882      delete [] candidateList;
883      delete [] givenDuals;
884      givenDuals=NULL;
885      candidateList=NULL;
886      int iRow;
887      for (iRow=0;iRow<numberRows_;iRow++) {
888        int iPivot=pivotVariable_[iRow];
889        clearPivoted(iPivot);
890      }
891    }
892  } else {
893    assert (!ifValuesPass);
894  }
895#ifdef CHECK_ACCURACY
896  {
897    if (numberIterations_) {
898      int il=-1;
899      double largest=1.0e-1;
900      int ilnb=-1;
901      double largestnb=1.0e-8;
902      for (int i=0;i<numberRows_+numberColumns_;i++) {
903        double diff = fabs(solution_[i]-zzzzzz[i]);
904        if (diff>largest) {
905          largest=diff;
906          il=i;
907        }
908        if (getColumnStatus(i)!=basic) {
909          if (diff>largestnb) {
910            largestnb=diff;
911            ilnb=i;
912          }
913        }
914      }
915      if (il>=0&&ilnb<0)
916        printf("largest diff of %g at %d, nonbasic %g at %d\n",
917               largest,il,largestnb,ilnb);
918    }
919  }
920#endif
921  while (problemStatus_==-1) {
922#ifdef CLP_DEBUG
923    if (givenDuals) {
924      double value5=0.0;
925      int i;
926      for (i=0;i<numberRows_+numberColumns_;i++) {
927        if (dj_[i]<-1.0e-6)
928          if (upper_[i]<1.0e20)
929            value5 += dj_[i]*upper_[i];
930          else
931            printf("bad dj %g on %d with large upper status %d\n",
932                   dj_[i],i,status_[i]&7);
933        else if (dj_[i] >1.0e-6)
934          if (lower_[i]>-1.0e20)
935            value5 += dj_[i]*lower_[i];
936          else
937            printf("bad dj %g on %d with large lower status %d\n",
938                   dj_[i],i,status_[i]&7);
939      }
940      printf("Values objective Value %g\n",value5);
941    }
942    if ((handler_->logLevel()&32)&&wasInValuesPass) {
943      double value5=0.0;
944      int i;
945      for (i=0;i<numberRows_+numberColumns_;i++) {
946        if (dj_[i]<-1.0e-6)
947          if (upper_[i]<1.0e20)
948            value5 += dj_[i]*upper_[i];
949        else if (dj_[i] >1.0e-6)
950          if (lower_[i]>-1.0e20)
951            value5 += dj_[i]*lower_[i];
952      }
953      printf("Values objective Value %g\n",value5);
954      {
955        int i;
956        for (i=0;i<numberRows_+numberColumns_;i++) {
957          int iSequence = i;
958          double oldValue;
959
960          switch(getStatus(iSequence)) {
961
962          case basic:
963          case ClpSimplex::isFixed:
964            break;
965          case isFree:
966          case superBasic:
967            abort();
968            break;
969          case atUpperBound:
970            oldValue = dj_[iSequence];
971            //assert (oldValue<=tolerance);
972            assert (fabs(solution_[iSequence]-upper_[iSequence])<1.0e-7);
973            break;
974          case atLowerBound:
975            oldValue = dj_[iSequence];
976            //assert (oldValue>=-tolerance);
977            assert (fabs(solution_[iSequence]-lower_[iSequence])<1.0e-7);
978            break;
979          }
980        }
981      }
982    }
983#endif
984#ifdef CLP_DEBUG
985    {
986      int i;
987      for (i=0;i<4;i++) {
988        rowArray_[i]->checkClear();
989      }
990      for (i=0;i<2;i++) {
991        columnArray_[i]->checkClear();
992      }
993    }
994#endif
995#if CLP_DEBUG>2
996    // very expensive
997    if (numberIterations_>3063&&numberIterations_<30700) {
998      //handler_->setLogLevel(63);
999      double saveValue = objectiveValue_;
1000      double * saveRow1 = new double[numberRows_];
1001      double * saveRow2 = new double[numberRows_];
1002      CoinMemcpyN(rowReducedCost_,numberRows_,saveRow1);
1003      CoinMemcpyN(rowActivityWork_,numberRows_,saveRow2);
1004      double * saveColumn1 = new double[numberColumns_];
1005      double * saveColumn2 = new double[numberColumns_];
1006      CoinMemcpyN(reducedCostWork_,numberColumns_,saveColumn1);
1007      CoinMemcpyN(columnActivityWork_,numberColumns_,saveColumn2);
1008      gutsOfSolution(NULL,NULL);
1009      printf("xxx %d old obj %g, recomputed %g, sum dual inf %g\n",
1010             numberIterations_,
1011             saveValue,objectiveValue_,sumDualInfeasibilities_);
1012      if (saveValue>objectiveValue_+1.0e-2)
1014      CoinMemcpyN(saveRow1,numberRows_,rowReducedCost_);
1015      CoinMemcpyN(saveRow2,numberRows_,rowActivityWork_);
1016      CoinMemcpyN(saveColumn1,numberColumns_,reducedCostWork_);
1017      CoinMemcpyN(saveColumn2,numberColumns_,columnActivityWork_);
1018      delete [] saveRow1;
1019      delete [] saveRow2;
1020      delete [] saveColumn1;
1021      delete [] saveColumn2;
1022      objectiveValue_=saveValue;
1023    }
1024#endif
1025#if 0
1026    //    if (factorization_->pivots()){
1027    {
1028      int iPivot;
1029      double * array = rowArray_[3]->denseVector();
1030      int i;
1031      for (iPivot=0;iPivot<numberRows_;iPivot++) {
1032        int iSequence = pivotVariable_[iPivot];
1033        unpack(rowArray_[3],iSequence);
1034        factorization_->updateColumn(rowArray_[2],rowArray_[3]);
1035        assert (fabs(array[iPivot]-1.0)<1.0e-4);
1036        array[iPivot]=0.0;
1037        for (i=0;i<numberRows_;i++)
1038          assert (fabs(array[i])<1.0e-4);
1039        rowArray_[3]->clear();
1040      }
1041    }
1042#endif
1043#ifdef CLP_DEBUG
1044    {
1045      int iSequence, number=numberRows_+numberColumns_;
1046      for (iSequence=0;iSequence<number;iSequence++) {
1047        double lowerValue=lower_[iSequence];
1048        double upperValue=upper_[iSequence];
1049        double value=solution_[iSequence];
1050        if(getStatus(iSequence)!=basic&&getStatus(iSequence)!=isFree) {
1051          assert(lowerValue>-1.0e20);
1052          assert(upperValue<1.0e20);
1053        }
1054        switch(getStatus(iSequence)) {
1055
1056        case basic:
1057          break;
1058        case isFree:
1059        case superBasic:
1060          break;
1061        case atUpperBound:
1062          assert (fabs(value-upperValue)<=primalTolerance_) ;
1063          break;
1064        case atLowerBound:
1065        case ClpSimplex::isFixed:
1066          assert (fabs(value-lowerValue)<=primalTolerance_) ;
1067          break;
1068        }
1069      }
1070    }
1071    if(numberIterations_==debugIteration) {
1072      printf("dodgy iteration coming up\n");
1073    }
1074#endif
1075    // choose row to go out
1076    // dualRow will go to virtual row pivot choice algorithm
1077    // make sure values pass off if it should be
1078    if (numberCandidates)
1079      candidate = candidateList[--numberCandidates];
1080    else
1081      candidate=-1;
1082    dualRow(candidate);
1083    if (pivotRow_>=0) {
1084      // we found a pivot row
1085      if (handler_->detail(CLP_SIMPLEX_PIVOTROW,messages_)<100) {
1086        handler_->message(CLP_SIMPLEX_PIVOTROW,messages_)
1087          <<pivotRow_
1088          <<CoinMessageEol;
1089      }
1090      // check accuracy of weights
1091      dualRowPivot_->checkAccuracy();
1092      // Get good size for pivot
1093      // Allow first few iterations to take tiny
1094      double acceptablePivot=1.0e-1*acceptablePivot_;
1095      if (numberIterations_>100)
1096        acceptablePivot=acceptablePivot_;
1097      if (factorization_->pivots()>10||
1098          (factorization_->pivots()&&saveSumDual))
1099        acceptablePivot=1.0e+3*acceptablePivot_; // if we have iterated be more strict
1100      else if (factorization_->pivots()>5)
1101        acceptablePivot=1.0e+2*acceptablePivot_; // if we have iterated be slightly more strict
1102      else if (factorization_->pivots())
1103        acceptablePivot=acceptablePivot_; // relax
1104      double bestPossiblePivot=1.0;
1105      // get sign for finding row of tableau
1106      if (candidate<0) {
1107        // normal iteration
1108        // create as packed
1109        double direction=directionOut_;
1110        rowArray_[0]->createPacked(1,&pivotRow_,&direction);
1111        factorization_->updateColumnTranspose(rowArray_[1],rowArray_[0]);
1112        // Allow to do dualColumn0
1114          spareIntArray_[0]=1;
1115        spareDoubleArray_[0]=acceptablePivot;
1116        rowArray_[3]->clear();
1117        sequenceIn_=-1;
1118        // put row of tableau in rowArray[0] and columnArray[0]
1119        if (!scaledMatrix_) {
1120          matrix_->transposeTimes(this,-1.0,
1121                                  rowArray_[0],rowArray_[3],columnArray_[0]);
1122        } else {
1123          double * saveR = rowScale_;
1124          double * saveC = columnScale_;
1125          rowScale_=NULL;
1126          columnScale_=NULL;
1127          scaledMatrix_->transposeTimes(this,-1.0,
1128                                  rowArray_[0],rowArray_[3],columnArray_[0]);
1129          rowScale_=saveR;
1130          columnScale_=saveC;
1131        }
1132        // do ratio test for normal iteration
1133        bestPossiblePivot = dualColumn(rowArray_[0],columnArray_[0],rowArray_[3],
1134                 columnArray_[1],acceptablePivot,dubiousWeights);
1135      } else {
1136        // Make sure direction plausible
1137        CoinAssert (upperOut_<1.0e50||lowerOut_>-1.0e50);
1138        // If in integer cleanup do direction using duals
1139        // may be wrong way round
1140        if(ifValuesPass==2) {
1141          if (dual_[pivotRow_]>0.0) {
1142            // this will give a -1 in pivot row (as slacks are -1.0)
1143            directionOut_ = 1;
1144          } else {
1145            directionOut_ = -1;
1146          }
1147        }
1148        if (directionOut_<0&&fabs(valueOut_-upperOut_)>dualBound_+primalTolerance_) {
1149          if (fabs(valueOut_-upperOut_)>fabs(valueOut_-lowerOut_))
1150            directionOut_=1;
1151        } else if (directionOut_>0&&fabs(valueOut_-lowerOut_)>dualBound_+primalTolerance_) {
1152          if (fabs(valueOut_-upperOut_)<fabs(valueOut_-lowerOut_))
1153            directionOut_=-1;
1154        }
1155        double direction=directionOut_;
1156        rowArray_[0]->createPacked(1,&pivotRow_,&direction);
1157        factorization_->updateColumnTranspose(rowArray_[1],rowArray_[0]);
1158        // put row of tableau in rowArray[0] and columnArray[0]
1159        if (!scaledMatrix_) {
1160          matrix_->transposeTimes(this,-1.0,
1161                                  rowArray_[0],rowArray_[3],columnArray_[0]);
1162        } else {
1163          double * saveR = rowScale_;
1164          double * saveC = columnScale_;
1165          rowScale_=NULL;
1166          columnScale_=NULL;
1167          scaledMatrix_->transposeTimes(this,-1.0,
1168                                  rowArray_[0],rowArray_[3],columnArray_[0]);
1169          rowScale_=saveR;
1170          columnScale_=saveC;
1171        }
1172        acceptablePivot *= 10.0;
1173        // do ratio test
1174        if (ifValuesPass==1) {
1175          checkPossibleValuesMove(rowArray_[0],columnArray_[0],
1176                                  acceptablePivot);
1177        } else {
1178          checkPossibleCleanup(rowArray_[0],columnArray_[0],
1179                                  acceptablePivot);
1180          if (sequenceIn_<0) {
1181            rowArray_[0]->clear();
1182            columnArray_[0]->clear();
1183            continue; // can't do anything
1184          }
1185        }
1186
1187        // recompute true dualOut_
1188        if (directionOut_<0) {
1189          dualOut_ = valueOut_ - upperOut_;
1190        } else {
1191          dualOut_ = lowerOut_ - valueOut_;
1192        }
1193        // check what happened if was values pass
1194        // may want to move part way i.e. movement
1195        bool normalIteration = (sequenceIn_!=sequenceOut_);
1196
1197        clearPivoted(sequenceOut_);  // make sure won't be done again
1198        // see if end of values pass
1199        if (!numberCandidates) {
1200          int iRow;
1201          delete [] candidateList;
1202          delete [] givenDuals;
1203          candidate=-2; // -2 signals end
1204          givenDuals=NULL;
1205          candidateList=NULL;
1206          ifValuesPass=1;
1207          for (iRow=0;iRow<numberRows_;iRow++) {
1208            int iPivot=pivotVariable_[iRow];
1209            //assert (fabs(dj_[iPivot]),1.0e-5);
1210            clearPivoted(iPivot);
1211          }
1212        }
1213        if (!normalIteration) {
1214          //rowArray_[0]->cleanAndPackSafe(1.0e-60);
1215          //columnArray_[0]->cleanAndPackSafe(1.0e-60);
1216          updateDualsInValuesPass(rowArray_[0],columnArray_[0],theta_);
1217          if (candidate==-2)
1218            problemStatus_=-2;
1219          continue; // skip rest of iteration
1220        } else {
1221          // recompute dualOut_
1222          if (directionOut_<0) {
1223            dualOut_ = valueOut_ - upperOut_;
1224          } else {
1225            dualOut_ = lowerOut_ - valueOut_;
1226          }
1227        }
1228      }
1229      if (sequenceIn_>=0) {
1230        // normal iteration
1231        // update the incoming column
1232        double btranAlpha = -alpha_*directionOut_; // for check
1233        unpackPacked(rowArray_[1]);
1234        // moved into updateWeights - factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
1235        // and update dual weights (can do in parallel - with extra array)
1236        alpha_ = dualRowPivot_->updateWeights(rowArray_[0],
1237                                              rowArray_[2],
1238                                              rowArray_[3],
1239                                              rowArray_[1]);
1240        // see if update stable
1241#ifdef CLP_DEBUG
1242        if ((handler_->logLevel()&32))
1243          printf("btran alpha %g, ftran alpha %g\n",btranAlpha,alpha_);
1244#endif
1245        double checkValue=1.0e-7;
1246        // if can't trust much and long way from optimal then relax
1247        if (largestPrimalError_>10.0)
1248          checkValue = CoinMin(1.0e-4,1.0e-8*largestPrimalError_);
1249        if (fabs(btranAlpha)<1.0e-12||fabs(alpha_)<1.0e-12||
1250            fabs(btranAlpha-alpha_)>checkValue*(1.0+fabs(alpha_))) {
1251          handler_->message(CLP_DUAL_CHECK,messages_)
1252            <<btranAlpha
1253            <<alpha_
1254            <<CoinMessageEol;
1255          if (factorization_->pivots()) {
1256            dualRowPivot_->unrollWeights();
1257            problemStatus_=-2; // factorize now
1258            rowArray_[0]->clear();
1259            rowArray_[1]->clear();
1260            columnArray_[0]->clear();
1261            returnCode=-2;
1262            break;
1263          } else {
1264            // take on more relaxed criterion
1265            double test;
1266            if (fabs(btranAlpha)<1.0e-8||fabs(alpha_)<1.0e-8)
1267              test = 1.0e-1*fabs(alpha_);
1268            else
1269              test = 1.0e-4*(1.0+fabs(alpha_));
1270            if (fabs(btranAlpha)<1.0e-12||fabs(alpha_)<1.0e-12||
1271                fabs(btranAlpha-alpha_)>test) {
1272              dualRowPivot_->unrollWeights();
1273              // need to reject something
1274              char x = isColumn(sequenceOut_) ? 'C' :'R';
1275              handler_->message(CLP_SIMPLEX_FLAG,messages_)
1276                <<x<<sequenceWithin(sequenceOut_)
1277                <<CoinMessageEol;
1278#ifdef COIN_DEVELOP
1279              printf("flag a %g %g\n",btranAlpha,alpha_);
1280#endif
1281              setFlagged(sequenceOut_);
1283              lastBadIteration_ = numberIterations_; // say be more cautious
1284              rowArray_[0]->clear();
1285              rowArray_[1]->clear();
1286              columnArray_[0]->clear();
1287              if (fabs(alpha_)<1.0e-10&&fabs(btranAlpha)<1.0e-8&&numberIterations_>100) {
1288                //printf("I think should declare infeasible\n");
1289                problemStatus_=1;
1290                returnCode=1;
1291                break;
1292              }
1293              continue;
1294            }
1295          }
1296        }
1297        // update duals BEFORE replaceColumn so can do updateColumn
1298        double objectiveChange=0.0;
1299        // do duals first as variables may flip bounds
1300        // rowArray_[0] and columnArray_[0] may have flips
1301        // so use rowArray_[3] for work array from here on
1302        int nswapped = 0;
1303        //rowArray_[0]->cleanAndPackSafe(1.0e-60);
1304        //columnArray_[0]->cleanAndPackSafe(1.0e-60);
1305        if (candidate==-1) {
1306          // make sure incoming doesn't count
1307          Status saveStatus = getStatus(sequenceIn_);
1308          setStatus(sequenceIn_,basic);
1309          nswapped = updateDualsInDual(rowArray_[0],columnArray_[0],
1310                                       rowArray_[2],theta_,
1311                                       objectiveChange,false);
1312          setStatus(sequenceIn_,saveStatus);
1313        } else {
1314          updateDualsInValuesPass(rowArray_[0],columnArray_[0],theta_);
1315        }
1316        double oldDualOut = dualOut_;
1317        // which will change basic solution
1318        if (nswapped) {
1319          factorization_->updateColumn(rowArray_[3],rowArray_[2]);
1320          dualRowPivot_->updatePrimalSolution(rowArray_[2],
1321                                              1.0,objectiveChange);
1322          // recompute dualOut_
1323          valueOut_ = solution_[sequenceOut_];
1324          if (directionOut_<0) {
1325            dualOut_ = valueOut_ - upperOut_;
1326          } else {
1327            dualOut_ = lowerOut_ - valueOut_;
1328          }
1329#if 0
1330          if (dualOut_<0.0) {
1331#ifdef CLP_DEBUG
1332            if (handler_->logLevel()&32) {
1333              printf(" dualOut_ %g %g save %g\n",dualOut_,averagePrimalInfeasibility,saveDualOut);
1334              printf("values %g %g %g %g %g %g %g\n",lowerOut_,valueOut_,upperOut_,
1335                     objectiveChange,);
1336            }
1337#endif
1338            if (upperOut_==lowerOut_)
1339              dualOut_=0.0;
1340          }
1341          if(dualOut_<-CoinMax(1.0e-12*averagePrimalInfeasibility,1.0e-8)
1342             &&factorization_->pivots()>100&&
1343             getStatus(sequenceIn_)!=isFree) {
1344            // going backwards - factorize
1345            dualRowPivot_->unrollWeights();
1346            problemStatus_=-2; // factorize now
1347            returnCode=-2;
1348            break;
1349          }
1350#endif
1351        }
1352        // amount primal will move
1353        double movement = -dualOut_*directionOut_/alpha_;
1354        double movementOld = oldDualOut*directionOut_/alpha_;
1355        // so objective should increase by fabs(dj)*movement
1356        // but we already have objective change - so check will be good
1357        if (objectiveChange+fabs(movementOld*dualIn_)<-CoinMax(1.0e-5,1.0e-12*fabs(objectiveValue_))) {
1358#ifdef CLP_DEBUG
1359          if (handler_->logLevel()&32)
1360            printf("movement %g, swap change %g, rest %g  * %g\n",
1361                   objectiveChange+fabs(movement*dualIn_),
1362                   objectiveChange,movement,dualIn_);
1363#endif
1364          if(factorization_->pivots()) {
1365            // going backwards - factorize
1366            dualRowPivot_->unrollWeights();
1367            problemStatus_=-2; // factorize now
1368            returnCode=-2;
1369            break;
1370          }
1371        }
1372        CoinAssert(fabs(dualOut_)<1.0e50);
1373        // if stable replace in basis
1375                                                         rowArray_[2],
1376                                                         rowArray_[1],
1377                                                         pivotRow_,
1378                                                         alpha_,
1379                                                     (moreSpecialOptions_&16)!=0);
1380        // if no pivots, bad update but reasonable alpha - take and invert
1382                   !factorization_->pivots()&&fabs(alpha_)>1.0e-5)
1385          // slight error
1387            problemStatus_=-2; // factorize now
1388            returnCode=-3;
1389          }
1390        } else if (updateStatus==2) {
1391          // major error
1392          dualRowPivot_->unrollWeights();
1393          // later we may need to unwind more e.g. fake bounds
1394          if (factorization_->pivots()&&
1395              ((moreSpecialOptions_&16)==0||factorization_->pivots()>4)) {
1396            problemStatus_=-2; // factorize now
1397            returnCode=-2;
1398            moreSpecialOptions_ |= 16;
1399            break;
1400          } else {
1401            // need to reject something
1402            char x = isColumn(sequenceOut_) ? 'C' :'R';
1403            handler_->message(CLP_SIMPLEX_FLAG,messages_)
1404              <<x<<sequenceWithin(sequenceOut_)
1405              <<CoinMessageEol;
1406#ifdef COIN_DEVELOP
1407            printf("flag b %g\n",alpha_);
1408#endif
1409            setFlagged(sequenceOut_);
1411            lastBadIteration_ = numberIterations_; // say be more cautious
1412            rowArray_[0]->clear();
1413            rowArray_[1]->clear();
1414            columnArray_[0]->clear();
1415            // make sure dual feasible
1416            // look at all rows and columns
1417            double objectiveChange=0.0;
1418            updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[1],
1419                              0.0,objectiveChange,true);
1420            continue;
1421          }
1422        } else if (updateStatus==3) {
1423          // out of memory
1424          // increase space if not many iterations
1425          if (factorization_->pivots()<
1426              0.5*factorization_->maximumPivots()&&
1427              factorization_->pivots()<200)
1428            factorization_->areaFactor(
1429                                       factorization_->areaFactor() * 1.1);
1430          problemStatus_=-2; // factorize now
1431        } else if (updateStatus==5) {
1432          problemStatus_=-2; // factorize now
1433        }
1434        // update primal solution
1435        if (theta_<0.0&&candidate==-1) {
1436#ifdef CLP_DEBUG
1437          if (handler_->logLevel()&32)
1438            printf("negative theta %g\n",theta_);
1439#endif
1440          theta_=0.0;
1441        }
1442        // do actual flips
1443        flipBounds(rowArray_[0],columnArray_[0],theta_);
1444        //rowArray_[1]->expand();
1445        dualRowPivot_->updatePrimalSolution(rowArray_[1],
1446                                            movement,
1447                                            objectiveChange);
1448#ifdef CLP_DEBUG
1449        double oldobj=objectiveValue_;
1450#endif
1451        // modify dualout
1452        dualOut_ /= alpha_;
1453        dualOut_ *= -directionOut_;
1454        //setStatus(sequenceIn_,basic);
1455        dj_[sequenceIn_]=0.0;
1456        double oldValue=valueIn_;
1457        if (directionIn_==-1) {
1458          // as if from upper bound
1459          valueIn_ = upperIn_+dualOut_;
1460        } else {
1461          // as if from lower bound
1462          valueIn_ = lowerIn_+dualOut_;
1463        }
1464        objectiveChange += cost_[sequenceIn_]*(valueIn_-oldValue);
1465        // outgoing
1466        // set dj to zero unless values pass
1467        if (directionOut_>0) {
1468          valueOut_ = lowerOut_;
1469          if (candidate==-1)
1470            dj_[sequenceOut_] = theta_;
1471        } else {
1472          valueOut_ = upperOut_;
1473          if (candidate==-1)
1474            dj_[sequenceOut_] = -theta_;
1475        }
1476        solution_[sequenceOut_]=valueOut_;
1477        int whatNext=housekeeping(objectiveChange);
1478#ifdef VUB
1479        {
1480          if ((sequenceIn_<numberColumns_&&vub[sequenceIn_]>=0)||toVub[sequenceIn_]>=0||
1481              (sequenceOut_<numberColumns_&&vub[sequenceOut_]>=0)||toVub[sequenceOut_]>=0) {
1482            int inSequence = sequenceIn_;
1483            int inVub = -1;
1484            if (sequenceIn_<numberColumns_)
1485              inVub = vub[sequenceIn_];
1486            int inBack = toVub[inSequence];
1487            int inSlack=-1;
1488            if (inSequence>=numberColumns_&&inBack>=0) {
1489              inSlack = inSequence-numberColumns_;
1490              inSequence = inBack;
1491              inBack = toVub[inSequence];
1492            }
1493            if (inVub>=0)
1494              printf("Vub %d in ",inSequence);
1495            if (inBack>=0&&inSlack<0)
1496              printf("%d (descendent of %d) in ",inSequence,inBack);
1497            if (inSlack>=0)
1498              printf("slack for row %d -> %d (descendent of %d) in ",inSlack,inSequence,inBack);
1499            int outSequence = sequenceOut_;
1500            int outVub = -1;
1501            if (sequenceOut_<numberColumns_)
1502              outVub = vub[sequenceOut_];
1503            int outBack = toVub[outSequence];
1504            int outSlack=-1;
1505            if (outSequence>=numberColumns_&&outBack>=0) {
1506              outSlack = outSequence-numberColumns_;
1507              outSequence = outBack;
1508              outBack = toVub[outSequence];
1509            }
1510            if (outVub>=0)
1511              printf("Vub %d out ",outSequence);
1512            if (outBack>=0&&outSlack<0)
1513              printf("%d (descendent of %d) out ",outSequence,outBack);
1514            if (outSlack>=0)
1515              printf("slack for row %d -> %d (descendent of %d) out ",outSlack,outSequence,outBack);
1516            printf("\n");
1517          }
1518        }
1519#endif
1520#if 0
1521        if (numberIterations_>206033)
1522          handler_->setLogLevel(63);
1523        if (numberIterations_>210567)
1524          exit(77);
1525#endif
1526        if (!givenDuals&&ifValuesPass&&ifValuesPass!=2) {
1527          handler_->message(CLP_END_VALUES_PASS,messages_)
1528            <<numberIterations_;
1529          whatNext=1;
1530        }
1531#ifdef CHECK_ACCURACY
1532        if (whatNext) {
1533          CoinMemcpyN(solution_,(numberRows_+numberColumns_),zzzzzz);
1534        }
1535#endif
1536        //if (numberIterations_==1890)
1537        //whatNext=1;
1538        //if (numberIterations_>2000)
1539        //exit(77);
1540        // and set bounds correctly
1541        originalBound(sequenceIn_);
1542        changeBound(sequenceOut_);
1543#ifdef CLP_DEBUG
1544        if (objectiveValue_<oldobj-1.0e-5&&(handler_->logLevel()&16))
1545          printf("obj backwards %g %g\n",objectiveValue_,oldobj);
1546#endif
1547#if 0
1548        {
1549          for (int i=0;i<numberRows_+numberColumns_;i++) {
1550            FakeBound bound = getFakeBound(i);
1551            if (bound==ClpSimplexDual::upperFake) {
1552              assert (upper_[i]<1.0e20);
1553            } else if (bound==ClpSimplexDual::lowerFake) {
1554              assert (lower_[i]>-1.0e20);
1555            } else if (bound==ClpSimplexDual::bothFake) {
1556              assert (upper_[i]<1.0e20);
1557              assert (lower_[i]>-1.0e20);
1558            }
1559          }
1560        }
1561#endif
1562        if (whatNext==1||candidate==-2) {
1563          problemStatus_ =-2; // refactorize
1564        } else if (whatNext==2) {
1565          // maximum iterations or equivalent
1566          problemStatus_= 3;
1567          returnCode=3;
1568          break;
1569        }
1570        // Check event
1571        {
1572          int status = eventHandler_->event(ClpEventHandler::endOfIteration);
1573          if (status>=0) {
1574            problemStatus_=5;
1575            secondaryStatus_=ClpEventHandler::endOfIteration;
1576            returnCode=4;
1577            break;
1578          }
1579        }
1580      } else {
1581#ifdef CLP_INVESTIGATE
1582        z_thinks=1;
1583#endif
1584        // no incoming column is valid
1585        pivotRow_=-1;
1586#ifdef CLP_DEBUG
1587        if (handler_->logLevel()&32)
1588          printf("** no column pivot\n");
1589#endif
1590        if (factorization_->pivots()<2&&acceptablePivot_<=1.0e-8) {
1591          //&&goodAccuracy()) {
1592          // If not in branch and bound etc save ray
1593          if ((specialOptions_&(1024|4096))==0) {
1594            // create ray anyway
1595            delete [] ray_;
1596            ray_ = new double [ numberRows_];
1597            rowArray_[0]->expand(); // in case packed
1598            CoinMemcpyN(rowArray_[0]->denseVector(),numberRows_,ray_);
1599          }
1600          // If we have just factorized and infeasibility reasonable say infeas
1601          if (((specialOptions_&4096)!=0||bestPossiblePivot<1.0e-11)&&dualBound_>1.0e8) {
1602            if (valueOut_>upperOut_+1.0e-4||valueOut_<lowerOut_-1.0e-4
1603                || (specialOptions_&64)==0) {
1604              // say infeasible
1605              problemStatus_=1;
1606              // unless primal feasible!!!!
1607              //printf("%d %g %d %g\n",numberPrimalInfeasibilities_,sumPrimalInfeasibilities_,
1608              //   numberDualInfeasibilities_,sumDualInfeasibilities_);
1609              //#define TEST_CLP_NODE
1610#ifndef TEST_CLP_NODE
1611              // Should be correct - but ...
1612              int numberFake = numberAtFakeBound();
1613              double sumPrimal =  (!numberFake) ? 2.0e5 : sumPrimalInfeasibilities_;
1614              if (sumPrimalInfeasibilities_<1.0e-3||sumDualInfeasibilities_>1.0e-5||
1615                  (sumPrimal<1.0e5&&(specialOptions_&1024)!=0&&factorization_->pivots())) {
1616                if (sumPrimal>50.0&&factorization_->pivots()>2) {
1617                  problemStatus_=-4;
1618#ifdef COIN_DEVELOP
1619                  printf("status to -4 at %d - primalinf %g pivots %d\n",
1620                         __LINE__,sumPrimalInfeasibilities_,
1621                         factorization_->pivots());
1622#endif
1623                } else {
1624                  problemStatus_=10;
1625#if COIN_DEVELOP>1
1626                  printf("returning at %d - primal %d %g - dual %d %g fake %d weight %g - pivs %d - options (1024-16384) %d %d %d %d %d\n",
1627                         __LINE__,numberPrimalInfeasibilities_,
1628                         sumPrimalInfeasibilities_,
1629                         numberDualInfeasibilities_,sumDualInfeasibilities_,
1630                         numberFake_,dualBound_,factorization_->pivots(),
1631                         (specialOptions_&1024)!=0 ? 1 : 0,
1632                         (specialOptions_&2048)!=0 ? 1 : 0,
1633                         (specialOptions_&4096)!=0 ? 1 : 0,
1634                         (specialOptions_&8192)!=0 ? 1 : 0,
1635                         (specialOptions_&16384)!=0 ? 1 : 0
1636                         );
1637#endif
1638                  // Get rid of objective
1639                  if ((specialOptions_&16384)==0)
1640                    objective_ = new ClpLinearObjective(NULL,numberColumns_);
1641                }
1642              }
1643#else
1644              if (sumPrimalInfeasibilities_<1.0e-3||sumDualInfeasibilities_>1.0e-6) {
1645#ifdef COIN_DEVELOP
1646                printf("at %d - primal %d %g - dual %d %g fake %d weight %g - pivs %d\n",
1647                       __LINE__,numberPrimalInfeasibilities_,
1648                       sumPrimalInfeasibilities_,
1649                       numberDualInfeasibilities_,sumDualInfeasibilities_,
1650                       numberFake_,dualBound_,factorization_->pivots());
1651#endif
1652                if ((specialOptions_&1024)!=0&&factorization_->pivots()) {
1653                  problemStatus_=10;
1654#if COIN_DEVELOP>1
1655                  printf("returning at %d\n",__LINE__);
1656#endif
1657                  // Get rid of objective
1658                  if ((specialOptions_&16384)==0)
1659                    objective_ = new ClpLinearObjective(NULL,numberColumns_);
1660                }
1661              }
1662#endif
1663              rowArray_[0]->clear();
1664              columnArray_[0]->clear();
1665              returnCode=1;
1666              break;
1667            }
1668          }
1669          // If special option set - put off as long as possible
1670          if ((specialOptions_&64)==0) {
1671            if (factorization_->pivots()==0)
1672              problemStatus_=-4; //say looks infeasible
1673          } else {
1674            // flag
1675            char x = isColumn(sequenceOut_) ? 'C' :'R';
1676            handler_->message(CLP_SIMPLEX_FLAG,messages_)
1677              <<x<<sequenceWithin(sequenceOut_)
1678              <<CoinMessageEol;
1679#ifdef COIN_DEVELOP
1680            printf("flag c\n");
1681#endif
1682            setFlagged(sequenceOut_);
1683            if (!factorization_->pivots()) {
1684              rowArray_[0]->clear();
1685              columnArray_[0]->clear();
1686              continue;
1687            }
1688          }
1689        }
1690        if (factorization_->pivots()<5&&acceptablePivot_>1.0e-8)
1691          acceptablePivot_=1.0e-8;
1692        rowArray_[0]->clear();
1693        columnArray_[0]->clear();
1694        returnCode=1;
1695        break;
1696      }
1697    } else {
1698#ifdef CLP_INVESTIGATE
1699      z_thinks=0;
1700#endif
1701      // no pivot row
1702#ifdef CLP_DEBUG
1703      if (handler_->logLevel()&32)
1704        printf("** no row pivot\n");
1705#endif
1706      // If in branch and bound try and get rid of fixed variables
1707      if ((specialOptions_&1024)!=0&&CLEAN_FIXED) {
1708        assert (!candidateList);
1709        candidateList = new int[numberRows_];
1710        int iRow;
1711        for (iRow=0;iRow<numberRows_;iRow++) {
1712          int iPivot=pivotVariable_[iRow];
1713          if (flagged(iPivot)||!pivoted(iPivot))
1714            continue;
1715          assert (iPivot<numberColumns_&&lower_[iPivot]==upper_[iPivot]);
1716          candidateList[numberCandidates++]=iRow;
1717        }
1718        // and set first candidate
1719        if (!numberCandidates) {
1720          delete [] candidateList;
1721          candidateList=NULL;
1722          int iRow;
1723          for (iRow=0;iRow<numberRows_;iRow++) {
1724            int iPivot=pivotVariable_[iRow];
1725            clearPivoted(iPivot);
1726          }
1727        } else {
1728          ifValuesPass=2;
1729          continue;
1730        }
1731      }
1732      int numberPivots = factorization_->pivots();
1733      bool specialCase;
1734      int useNumberFake;
1735      returnCode=0;
1736      if (numberPivots<=CoinMax(dontFactorizePivots_,20)&&
1737          (specialOptions_&2048)!=0&&(true||!numberChanged_||perturbation_==101)
1738          &&dualBound_>=1.0e8) {
1739        specialCase=true;
1740        // as dual bound high - should be okay
1741        useNumberFake=0;
1742      } else {
1743        specialCase=false;
1744        useNumberFake=numberFake_;
1745      }
1746      if (!numberPivots||specialCase) {
1747        // may have crept through - so may be optimal
1748        // check any flagged variables
1749        int iRow;
1750        for (iRow=0;iRow<numberRows_;iRow++) {
1751          int iPivot=pivotVariable_[iRow];
1752          if (flagged(iPivot))
1753            break;
1754        }
1755        if (iRow<numberRows_&&numberPivots) {
1756          // try factorization
1757          returnCode=-2;
1758        }
1759
1760        if (useNumberFake||numberDualInfeasibilities_) {
1761          // may be dual infeasible
1762          if ((specialOptions_&1024)==0)
1763            problemStatus_=-5;
1764        } else {
1765          if (iRow<numberRows_) {
1766#ifdef COIN_DEVELOP
1767            std::cout<<"Flagged variables at end - infeasible?"<<std::endl;
1768            printf("Probably infeasible - pivot was %g\n",alpha_);
1769#endif
1770            //if (fabs(alpha_)<1.0e-4) {
1771            //problemStatus_=1;
1772            //} else {
1773#ifdef CLP_DEBUG
1774            abort();
1775#endif
1776            //}
1777            problemStatus_=-5;
1778          } else {
1779            problemStatus_=0;
1780            // May be perturbed
1781            if (perturbation_==101||numberChanged_) {
1782              perturbation_=102; // stop any perturbations
1783              //double changeCost;
1784              //changeBounds(1,NULL,changeCost);
1785              createRim4(false);
1786              // make sure duals are current
1787              computeDuals(givenDuals);
1788              checkDualSolution();
1789              if (numberDualInfeasibilities_) {
1790                problemStatus_=10; // was -3;
1791              } else {
1792                computeObjectiveValue(true);
1793              }
1794            } else if (numberPivots) {
1795              computeObjectiveValue(true);
1796            }
1797            if (numberPivots<-1000) {
1798              // objective may be wrong
1799              objectiveValue_ = innerProduct(cost_,numberColumns_+numberRows_,solution_);
1800              objectiveValue_ += objective_->nonlinearOffset();
1801              objectiveValue_ /= (objectiveScale_*rhsScale_);
1802              if ((specialOptions_&16384)==0) {
1803                // and dual_ may be wrong (i.e. for fixed or basic)
1804                CoinIndexedVector * arrayVector = rowArray_[1];
1805                arrayVector->clear();
1806                int iRow;
1807                double * array = arrayVector->denseVector();
1808                /* Use dual_ instead of array
1809                   Even though dual_ is only numberRows_ long this is
1810                   okay as gets permuted to longer rowArray_[2]
1811                */
1812                arrayVector->setDenseVector(dual_);
1813                int * index = arrayVector->getIndices();
1814                int number=0;
1815                for (iRow=0;iRow<numberRows_;iRow++) {
1816                  int iPivot=pivotVariable_[iRow];
1817                  double value = cost_[iPivot];
1818                  dual_[iRow]=value;
1819                  if (value) {
1820                    index[number++]=iRow;
1821                  }
1822                }
1823                arrayVector->setNumElements(number);
1824                // Extended duals before "updateTranspose"
1825                matrix_->dualExpanded(this,arrayVector,NULL,0);
1826                // Btran basic costs
1827                rowArray_[2]->clear();
1828                factorization_->updateColumnTranspose(rowArray_[2],arrayVector);
1829                // and return vector
1830                arrayVector->setDenseVector(array);
1831              }
1832            }
1833            sumPrimalInfeasibilities_=0.0;
1834            if ((specialOptions_&(1024+16384))!=0&&!problemStatus_) {
1835              CoinIndexedVector * arrayVector = rowArray_[1];
1836              arrayVector->clear();
1837              double * rhs = arrayVector->denseVector();
1838              times(1.0,solution_,rhs);
1839#ifdef CHECK_ACCURACY
1841#endif
1843              int i;
1844              for ( i=0;i<numberRows_;i++) {
1845                if (rhs[i]<rowLowerWork_[i]-primalTolerance_||
1846                    rhs[i]>rowUpperWork_[i]+primalTolerance_) {
1848#ifdef CHECK_ACCURACY
1849                  printf("row %d out of bounds %g, %g correct %g bad %g\n",i,
1850                         rowLowerWork_[i],rowUpperWork_[i],
1851                         rhs[i],rowActivityWork_[i]);
1852#endif
1853                } else if (fabs(rhs[i]-rowActivityWork_[i])>1.0e-3) {
1854#ifdef CHECK_ACCURACY
1856                  printf("row %d correct %g bad %g\n",i,rhs[i],rowActivityWork_[i]);
1857#endif
1858                }
1859                rhs[i]=0.0;
1860              }
1861              for ( i=0;i<numberColumns_;i++) {
1862                if (solution_[i]<columnLowerWork_[i]-primalTolerance_||
1863                    solution_[i]>columnUpperWork_[i]+primalTolerance_) {
1865#ifdef CHECK_ACCURACY
1866                  printf("column %d out of bounds %g, %g correct %g bad %g\n",i,
1867                         columnLowerWork_[i],columnUpperWork_[i],
1868                         solution_[i],columnActivityWork_[i]);
1869#endif
1870                }
1871              }
1873                problemStatus_=-3;
1874                returnCode=-2;
1875                // Force to re-factorize early next time
1876                int numberPivots = factorization_->pivots();
1877                forceFactorization_=CoinMin(forceFactorization_,(numberPivots+1)>>1);
1878              }
1879            }
1880          }
1881        }
1882      } else {
1883        problemStatus_=-3;
1884        returnCode=-2;
1885        // Force to re-factorize early next time
1886        int numberPivots = factorization_->pivots();
1887        forceFactorization_=CoinMin(forceFactorization_,(numberPivots+1)>>1);
1888      }
1889      break;
1890    }
1891  }
1892  if (givenDuals) {
1893    CoinMemcpyN(dj_,numberRows_+numberColumns_,givenDuals);
1894    // get rid of any values pass array
1895    delete [] candidateList;
1896  }
1897  delete [] dubiousWeights;
1898  return returnCode;
1899}
1900/* The duals are updated by the given arrays.
1901   Returns number of infeasibilities.
1902   rowArray and columnarray will have flipped
1903   The output vector has movement (row length array) */
1904int
1905ClpSimplexDual::updateDualsInDual(CoinIndexedVector * rowArray,
1906                                  CoinIndexedVector * columnArray,
1907                                  CoinIndexedVector * outputArray,
1908                                  double theta,
1909                                  double & objectiveChange,
1910                                  bool fullRecompute)
1911{
1912
1913  outputArray->clear();
1914
1915
1916  int numberInfeasibilities=0;
1917  int numberRowInfeasibilities=0;
1918
1919  // get a tolerance
1920  double tolerance = dualTolerance_;
1921  // we can't really trust infeasibilities if there is dual error
1922  double error = CoinMin(1.0e-2,largestDualError_);
1923  // allow tolerance at least slightly bigger than standard
1924  tolerance = tolerance +  error;
1925
1926  double changeObj=0.0;
1927
1928  // Coding is very similar but we can save a bit by splitting
1929  // Do rows
1930  if (!fullRecompute) {
1931    int i;
1932    double * reducedCost = djRegion(0);
1933    const double * lower = lowerRegion(0);
1934    const double * upper = upperRegion(0);
1935    const double * cost = costRegion(0);
1936    double * work;
1937    int number;
1938    int * which;
1939    assert(rowArray->packedMode());
1940    work = rowArray->denseVector();
1941    number = rowArray->getNumElements();
1942    which = rowArray->getIndices();
1943    for (i=0;i<number;i++) {
1944      int iSequence = which[i];
1945      double alphaI = work[i];
1946      work[i]=0.0;
1947
1948      Status status = getStatus(iSequence+numberColumns_);
1949      // more likely to be at upper bound ?
1950      if (status==atUpperBound) {
1951
1952        double value = reducedCost[iSequence]-theta*alphaI;
1953        reducedCost[iSequence]=value;
1954
1955        if (value>tolerance) {
1956          // to lower bound (if swap)
1957          which[numberInfeasibilities++]=iSequence;
1958          double movement = lower[iSequence]-upper[iSequence];
1959          assert (fabs(movement)<1.0e30);
1960#ifdef CLP_DEBUG
1961          if ((handler_->logLevel()&32))
1962            printf("%d %d, new dj %g, alpha %g, movement %g\n",
1963                   0,iSequence,value,alphaI,movement);
1964#endif
1965          changeObj += movement*cost[iSequence];
1967        }
1968      } else if (status==atLowerBound) {
1969
1970        double value = reducedCost[iSequence]-theta*alphaI;
1971        reducedCost[iSequence]=value;
1972
1973        if (value<-tolerance) {
1974          // to upper bound
1975          which[numberInfeasibilities++]=iSequence;
1976          double movement = upper[iSequence] - lower[iSequence];
1977          assert (fabs(movement)<1.0e30);
1978#ifdef CLP_DEBUG
1979          if ((handler_->logLevel()&32))
1980            printf("%d %d, new dj %g, alpha %g, movement %g\n",
1981                   0,iSequence,value,alphaI,movement);
1982#endif
1983          changeObj += movement*cost[iSequence];
1985        }
1986      }
1987    }
1988  } else  {
1989    double * solution = solutionRegion(0);
1990    double * reducedCost = djRegion(0);
1991    const double * lower = lowerRegion(0);
1992    const double * upper = upperRegion(0);
1993    const double * cost = costRegion(0);
1994    int * which;
1995    which = rowArray->getIndices();
1996    int iSequence;
1997    for (iSequence=0;iSequence<numberRows_;iSequence++) {
1998      double value = reducedCost[iSequence];
1999
2000      Status status = getStatus(iSequence+numberColumns_);
2001      // more likely to be at upper bound ?
2002      if (status==atUpperBound) {
2003        double movement=0.0;
2004
2005        if (value>tolerance) {
2006          // to lower bound (if swap)
2007          // put back alpha
2008          which[numberInfeasibilities++]=iSequence;
2009          movement = lower[iSequence]-upper[iSequence];
2010          changeObj += movement*cost[iSequence];
2012          assert (fabs(movement)<1.0e30);
2013        } else if (value>-tolerance) {
2014          // at correct bound but may swap
2015          FakeBound bound = getFakeBound(iSequence+numberColumns_);
2016          if (bound==ClpSimplexDual::upperFake) {
2017            movement = lower[iSequence]-upper[iSequence];
2018            assert (fabs(movement)<1.0e30);
2019            setStatus(iSequence+numberColumns_,atLowerBound);
2020            solution[iSequence] = lower[iSequence];
2021            changeObj += movement*cost[iSequence];
2022            numberFake_--;
2023            setFakeBound(iSequence+numberColumns_,noFake);
2024          }
2025        }
2026      } else if (status==atLowerBound) {
2027        double movement=0.0;
2028
2029        if (value<-tolerance) {
2030          // to upper bound
2031          // put back alpha
2032          which[numberInfeasibilities++]=iSequence;
2033          movement = upper[iSequence] - lower[iSequence];
2034          assert (fabs(movement)<1.0e30);
2035          changeObj += movement*cost[iSequence];
2037        } else if (value<tolerance) {
2038          // at correct bound but may swap
2039          FakeBound bound = getFakeBound(iSequence+numberColumns_);
2040          if (bound==ClpSimplexDual::lowerFake) {
2041            movement = upper[iSequence]-lower[iSequence];
2042            assert (fabs(movement)<1.0e30);
2043            setStatus(iSequence+numberColumns_,atUpperBound);
2044            solution[iSequence] = upper[iSequence];
2045            changeObj += movement*cost[iSequence];
2046            numberFake_--;
2047            setFakeBound(iSequence+numberColumns_,noFake);
2048          }
2049        }
2050      }
2051    }
2052  }
2053
2054  // Do columns
2055  if (!fullRecompute) {
2056    int i;
2057    double * reducedCost = djRegion(1);
2058    const double * lower = lowerRegion(1);
2059    const double * upper = upperRegion(1);
2060    const double * cost = costRegion(1);
2061    double * work;
2062    int number;
2063    int * which;
2064    // set number of infeasibilities in row array
2065    numberRowInfeasibilities=numberInfeasibilities;
2066    rowArray->setNumElements(numberInfeasibilities);
2067    numberInfeasibilities=0;
2068    work = columnArray->denseVector();
2069    number = columnArray->getNumElements();
2070    which = columnArray->getIndices();
2071
2072    for (i=0;i<number;i++) {
2073      int iSequence = which[i];
2074      double alphaI = work[i];
2075      work[i]=0.0;
2076
2077      Status status = getStatus(iSequence);
2078      if (status==atLowerBound) {
2079        double value = reducedCost[iSequence]-theta*alphaI;
2080        reducedCost[iSequence]=value;
2081        double movement=0.0;
2082
2083        if (value<-tolerance) {
2084          // to upper bound
2085          which[numberInfeasibilities++]=iSequence;
2086          movement = upper[iSequence] - lower[iSequence];
2087          assert (fabs(movement)<1.0e30);
2088#ifdef CLP_DEBUG
2089          if ((handler_->logLevel()&32))
2090            printf("%d %d, new dj %g, alpha %g, movement %g\n",
2091                   1,iSequence,value,alphaI,movement);
2092#endif
2093          changeObj += movement*cost[iSequence];
2095        }
2096      } else if (status==atUpperBound) {
2097        double value = reducedCost[iSequence]-theta*alphaI;
2098        reducedCost[iSequence]=value;
2099        double movement=0.0;
2100
2101        if (value>tolerance) {
2102          // to lower bound (if swap)
2103          which[numberInfeasibilities++]=iSequence;
2104          movement = lower[iSequence]-upper[iSequence];
2105          assert (fabs(movement)<1.0e30);
2106#ifdef CLP_DEBUG
2107          if ((handler_->logLevel()&32))
2108            printf("%d %d, new dj %g, alpha %g, movement %g\n",
2109                   1,iSequence,value,alphaI,movement);
2110#endif
2111          changeObj += movement*cost[iSequence];
2113        }
2114      } else if (status==isFree) {
2115        double value = reducedCost[iSequence]-theta*alphaI;
2116        reducedCost[iSequence]=value;
2117      }
2118    }
2119  } else {
2120    double * solution = solutionRegion(1);
2121    double * reducedCost = djRegion(1);
2122    const double * lower = lowerRegion(1);
2123    const double * upper = upperRegion(1);
2124    const double * cost = costRegion(1);
2125    int * which;
2126    // set number of infeasibilities in row array
2127    numberRowInfeasibilities=numberInfeasibilities;
2128    rowArray->setNumElements(numberInfeasibilities);
2129    numberInfeasibilities=0;
2130    which = columnArray->getIndices();
2131    int iSequence;
2132    for (iSequence=0;iSequence<numberColumns_;iSequence++) {
2133      double value = reducedCost[iSequence];
2134
2135      Status status = getStatus(iSequence);
2136      if (status==atLowerBound) {
2137        double movement=0.0;
2138
2139        if (value<-tolerance) {
2140          // to upper bound
2141          // put back alpha
2142          which[numberInfeasibilities++]=iSequence;
2143          movement = upper[iSequence] - lower[iSequence];
2144          assert (fabs(movement)<1.0e30);
2145          changeObj += movement*cost[iSequence];
2147        } else if (value<tolerance) {
2148          // at correct bound but may swap
2149          FakeBound bound = getFakeBound(iSequence);
2150          if (bound==ClpSimplexDual::lowerFake) {
2151            movement = upper[iSequence]-lower[iSequence];
2152            assert (fabs(movement)<1.0e30);
2153            setStatus(iSequence,atUpperBound);
2154            solution[iSequence] = upper[iSequence];
2155            changeObj += movement*cost[iSequence];
2156            numberFake_--;
2157            setFakeBound(iSequence,noFake);
2158          }
2159        }
2160      } else if (status==atUpperBound) {
2161        double movement=0.0;
2162
2163        if (value>tolerance) {
2164          // to lower bound (if swap)
2165          // put back alpha
2166          which[numberInfeasibilities++]=iSequence;
2167          movement = lower[iSequence]-upper[iSequence];
2168          assert (fabs(movement)<1.0e30);
2169          changeObj += movement*cost[iSequence];
2171        } else if (value>-tolerance) {
2172          // at correct bound but may swap
2173          FakeBound bound = getFakeBound(iSequence);
2174          if (bound==ClpSimplexDual::upperFake) {
2175            movement = lower[iSequence]-upper[iSequence];
2176            assert (fabs(movement)<1.0e30);
2177            setStatus(iSequence,atLowerBound);
2178            solution[iSequence] = lower[iSequence];
2179            changeObj += movement*cost[iSequence];
2180            numberFake_--;
2181            setFakeBound(iSequence,noFake);
2182          }
2183        }
2184      }
2185    }
2186  }
2187#ifdef CLP_DEBUG
2188  if (fullRecompute&&numberFake_&&(handler_->logLevel()&16)!=0)
2189    printf("%d fake after full update\n",numberFake_);
2190#endif
2191  // set number of infeasibilities
2192  columnArray->setNumElements(numberInfeasibilities);
2193  numberInfeasibilities += numberRowInfeasibilities;
2194  if (fullRecompute) {
2195    // do actual flips
2196    flipBounds(rowArray,columnArray,0.0);
2197  }
2198  objectiveChange += changeObj;
2199  return numberInfeasibilities;
2200}
2201void
2202ClpSimplexDual::updateDualsInValuesPass(CoinIndexedVector * rowArray,
2203                                  CoinIndexedVector * columnArray,
2204                                        double theta)
2205{
2206
2207  // use a tighter tolerance except for all being okay
2208  double tolerance = dualTolerance_;
2209
2210  // Coding is very similar but we can save a bit by splitting
2211  // Do rows
2212  {
2213    int i;
2214    double * reducedCost = djRegion(0);
2215    double * work;
2216    int number;
2217    int * which;
2218    work = rowArray->denseVector();
2219    number = rowArray->getNumElements();
2220    which = rowArray->getIndices();
2221    for (i=0;i<number;i++) {
2222      int iSequence = which[i];
2223      double alphaI = work[i];
2224      double value = reducedCost[iSequence]-theta*alphaI;
2225      work[i]=0.0;
2226      reducedCost[iSequence]=value;
2227
2228      Status status = getStatus(iSequence+numberColumns_);
2229      // more likely to be at upper bound ?
2230      if (status==atUpperBound) {
2231
2232        if (value>tolerance)
2233          reducedCost[iSequence]=0.0;
2234      } else if (status==atLowerBound) {
2235
2236        if (value<-tolerance) {
2237          reducedCost[iSequence]=0.0;
2238        }
2239      }
2240    }
2241  }
2242  rowArray->setNumElements(0);
2243
2244  // Do columns
2245  {
2246    int i;
2247    double * reducedCost = djRegion(1);
2248    double * work;
2249    int number;
2250    int * which;
2251    work = columnArray->denseVector();
2252    number = columnArray->getNumElements();
2253    which = columnArray->getIndices();
2254
2255    for (i=0;i<number;i++) {
2256      int iSequence = which[i];
2257      double alphaI = work[i];
2258      double value = reducedCost[iSequence]-theta*alphaI;
2259      work[i]=0.0;
2260      reducedCost[iSequence]=value;
2261
2262      Status status = getStatus(iSequence);
2263      if (status==atLowerBound) {
2264        if (value<-tolerance)
2265          reducedCost[iSequence]=0.0;
2266      } else if (status==atUpperBound) {
2267        if (value>tolerance)
2268          reducedCost[iSequence]=0.0;
2269      }
2270    }
2271  }
2272  columnArray->setNumElements(0);
2273}
2274/*
2275   Chooses dual pivot row
2276   Would be faster with separate region to scan
2277   and will have this (with square of infeasibility) when steepest
2278   For easy problems we can just choose one of the first rows we look at
2279*/
2280void
2282{
2283  // get pivot row using whichever method it is
2284  int chosenRow=-1;
2285#ifdef FORCE_FOLLOW
2286  bool forceThis=false;
2287  if (!fpFollow&&strlen(forceFile)) {
2289    assert (fpFollow);
2290  }
2291  if (fpFollow) {
2292    if (numberIterations_<=force_iteration) {
2293      // read to next Clp0102
2294      char temp[300];
2295      while (fgets(temp,250,fpFollow)) {
2296        if (strncmp(temp,"Clp0102",7))
2297          continue;
2298        char cin,cout;
2299        sscanf(temp+9,"%d%*f%*s%*c%c%d%*s%*c%c%d",
2300               &force_iteration,&cin,&force_in,&cout,&force_out);
2301        if (cin=='R')
2302          force_in += numberColumns_;
2303        if (cout=='R')
2304          force_out += numberColumns_;
2305        forceThis=true;
2306        assert (numberIterations_==force_iteration-1);
2307        printf("Iteration %d will force %d out and %d in\n",
2308               force_iteration,force_out,force_in);
2310        break;
2311      }
2312    } else {
2313      // use old
2314      forceThis=true;
2315    }
2316    if (!forceThis) {
2317      fclose(fpFollow);
2318      fpFollow=NULL;
2319      forceFile="";
2320    }
2321  }
2322#endif
2323  double freeAlpha=0.0;
2325    // first see if any free variables and put them in basis
2326    int nextFree = nextSuperBasic();
2327    //nextFree=-1; //off
2328    if (nextFree>=0) {
2329      // unpack vector and find a good pivot
2330      unpack(rowArray_[1],nextFree);
2331      factorization_->updateColumn(rowArray_[2],rowArray_[1]);
2332
2333      double * work=rowArray_[1]->denseVector();
2334      int number=rowArray_[1]->getNumElements();
2335      int * which=rowArray_[1]->getIndices();
2336      double bestFeasibleAlpha=0.0;
2337      int bestFeasibleRow=-1;
2338      double bestInfeasibleAlpha=0.0;
2339      int bestInfeasibleRow=-1;
2340      int i;
2341
2342      for (i=0;i<number;i++) {
2343        int iRow = which[i];
2344        double alpha = fabs(work[iRow]);
2345        if (alpha>1.0e-3) {
2346          int iSequence=pivotVariable_[iRow];
2347          double value = solution_[iSequence];
2348          double lower = lower_[iSequence];
2349          double upper = upper_[iSequence];
2350          double infeasibility=0.0;
2351          if (value>upper)
2352            infeasibility = value-upper;
2353          else if (value<lower)
2354            infeasibility = lower-value;
2355          if (infeasibility*alpha>bestInfeasibleAlpha&&alpha>1.0e-1) {
2356            if (!flagged(iSequence)) {
2357              bestInfeasibleAlpha = infeasibility*alpha;
2358              bestInfeasibleRow=iRow;
2359            }
2360          }
2361          if (alpha>bestFeasibleAlpha&&(lower>-1.0e20||upper<1.0e20)) {
2362            bestFeasibleAlpha = alpha;
2363            bestFeasibleRow=iRow;
2364          }
2365        }
2366      }
2367      if (bestInfeasibleRow>=0)
2368        chosenRow = bestInfeasibleRow;
2369      else if (bestFeasibleAlpha>1.0e-2)
2370        chosenRow = bestFeasibleRow;
2371      if (chosenRow>=0) {
2372        pivotRow_=chosenRow;
2373        freeAlpha=work[chosenRow];
2374      }
2375      rowArray_[1]->clear();
2376    }
2377  } else {
2378    // in values pass
2380#ifdef FORCE_FOLLOW
2381    if(forceThis) {
2383      chosenRow=-1;
2384      for (int i=0;i<numberRows_;i++) {
2385        if (pivotVariable_[i]==force_out) {
2386          chosenRow=i;
2387          break;
2388        }
2389      }
2390      assert (chosenRow>=0);
2391    }
2392#endif
2393    pivotRow_=chosenRow;
2394  }
2395  if (chosenRow<0)
2396    pivotRow_=dualRowPivot_->pivotRow();
2397
2398  if (pivotRow_>=0) {
2399    sequenceOut_ = pivotVariable_[pivotRow_];
2400    valueOut_ = solution_[sequenceOut_];
2401    lowerOut_ = lower_[sequenceOut_];
2402    upperOut_ = upper_[sequenceOut_];
2404      // if we have problems we could try other way and hope we get a
2405      // zero pivot?
2406      if (valueOut_>upperOut_) {
2407        directionOut_ = -1;
2408        dualOut_ = valueOut_ - upperOut_;
2409      } else if (valueOut_<lowerOut_) {
2410        directionOut_ = 1;
2411        dualOut_ = lowerOut_ - valueOut_;
2412      } else {
2413#if 1
2414        // odd (could be free) - it's feasible - go to nearest
2415        if (valueOut_-lowerOut_<upperOut_-valueOut_) {
2416          directionOut_ = 1;
2417          dualOut_ = lowerOut_ - valueOut_;
2418        } else {
2419          directionOut_ = -1;
2420          dualOut_ = valueOut_ - upperOut_;
2421        }
2422#else
2423        // odd (could be free) - it's feasible - improve obj
2424        printf("direction from alpha of %g is %d\n",
2425               freeAlpha,freeAlpha>0.0 ? 1 : -1);
2426        if (valueOut_-lowerOut_>1.0e20)
2427          freeAlpha=1.0;
2428        else if(upperOut_-valueOut_>1.0e20)
2429          freeAlpha=-1.0;
2430        //if (valueOut_-lowerOut_<upperOut_-valueOut_) {
2431        if (freeAlpha<0.0) {
2432          directionOut_ = 1;
2433          dualOut_ = lowerOut_ - valueOut_;
2434        } else {
2435          directionOut_ = -1;
2436          dualOut_ = valueOut_ - upperOut_;
2437        }
2438        printf("direction taken %d - bounds %g %g %g\n",
2439               directionOut_,lowerOut_,valueOut_,upperOut_);
2440#endif
2441      }
2442#ifdef CLP_DEBUG
2443      assert(dualOut_>=0.0);
2444#endif
2445    } else {
2446      // in values pass so just use sign of dj
2447      // We don't want to go through any barriers so set dualOut low
2448      // free variables will never be here
2449      dualOut_ = 1.0e-6;
2450      if (dj_[sequenceOut_]>0.0) {
2451        // this will give a -1 in pivot row (as slacks are -1.0)
2452        directionOut_ = 1;
2453      } else {
2454        directionOut_ = -1;
2455      }
2456    }
2457  }
2458  return ;
2459}
2460// Checks if any fake bounds active - if so returns number and modifies
2461// dualBound_ and everything.
2462// Free variables will be left as free
2463// Returns number of bounds changed if >=0
2464// Returns -1 if not initialize and no effect
2465// Fills in changeVector which can be used to see if unbounded
2466// and cost of change vector
2467int
2468ClpSimplexDual::changeBounds(int initialize,
2469                                 CoinIndexedVector * outputArray,
2470                                 double & changeCost)
2471{
2472  numberFake_ = 0;
2473  if (!initialize) {
2474    int numberInfeasibilities;
2475    double newBound;
2476    newBound = 5.0*dualBound_;
2477    numberInfeasibilities=0;
2478    changeCost=0.0;
2479    // put back original bounds and then check
2480    createRim1(false);
2481    int iSequence;
2482    // bounds will get bigger - just look at ones at bounds
2483    for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
2484      double lowerValue=lower_[iSequence];
2485      double upperValue=upper_[iSequence];
2486      double value=solution_[iSequence];
2487      setFakeBound(iSequence,ClpSimplexDual::noFake);
2488      switch(getStatus(iSequence)) {
2489
2490      case basic:
2491      case ClpSimplex::isFixed:
2492        break;
2493      case isFree:
2494      case superBasic:
2495        break;
2496      case atUpperBound:
2497        if (fabs(value-upperValue)>primalTolerance_)
2498          numberInfeasibilities++;
2499        break;
2500      case atLowerBound:
2501        if (fabs(value-lowerValue)>primalTolerance_)
2502          numberInfeasibilities++;
2503        break;
2504      }
2505    }
2506    // If dual infeasible then carry on
2507    if (numberInfeasibilities) {
2508      handler_->message(CLP_DUAL_CHECKB,messages_)
2509        <<newBound
2510        <<CoinMessageEol;
2511      int iSequence;
2512      for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
2513        double lowerValue=lower_[iSequence];
2514        double upperValue=upper_[iSequence];
2515        double newLowerValue;
2516        double newUpperValue;
2517        Status status = getStatus(iSequence);
2518        if (status==atUpperBound||
2519            status==atLowerBound) {
2520          double value = solution_[iSequence];
2521          if (value-lowerValue<=upperValue-value) {
2522            newLowerValue = CoinMax(lowerValue,value-0.666667*newBound);
2523            newUpperValue = CoinMin(upperValue,newLowerValue+newBound);
2524          } else {
2525            newUpperValue = CoinMin(upperValue,value+0.666667*newBound);
2526            newLowerValue = CoinMax(lowerValue,newUpperValue-newBound);
2527          }
2528          lower_[iSequence]=newLowerValue;
2529          upper_[iSequence]=newUpperValue;
2530          if (newLowerValue > lowerValue) {
2531            if (newUpperValue < upperValue) {
2532              setFakeBound(iSequence,ClpSimplexDual::bothFake);
2533              numberFake_++;
2534            } else {
2535              setFakeBound(iSequence,ClpSimplexDual::lowerFake);
2536              numberFake_++;
2537            }
2538          } else {
2539            if (newUpperValue < upperValue) {
2540              setFakeBound(iSequence,ClpSimplexDual::upperFake);
2541              numberFake_++;
2542            }
2543          }
2544          if (status==atUpperBound)
2545            solution_[iSequence] = newUpperValue;
2546          else
2547            solution_[iSequence] = newLowerValue;
2548          double movement = solution_[iSequence] - value;
2549          if (movement&&outputArray) {
2550            if (iSequence>=numberColumns_) {
2552              changeCost += movement*cost_[iSequence];
2553            } else {
2555              changeCost += movement*cost_[iSequence];
2556            }
2557          }
2558        }
2559      }
2560      dualBound_ = newBound;
2561    } else {
2562      numberInfeasibilities=-1;
2563    }
2564    return numberInfeasibilities;
2565  } else if (initialize==1) {
2566    int iSequence;
2567
2568    for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
2569      Status status = getStatus(iSequence);
2570      if (status==atUpperBound||
2571          status==atLowerBound) {
2572        double lowerValue=lower_[iSequence];
2573        double upperValue=upper_[iSequence];
2574        double value = solution_[iSequence];
2575        if (lowerValue>-largeValue_||upperValue<largeValue_) {
2576          if (lowerValue-value>-0.5*dualBound_||
2577              upperValue-value<0.5*dualBound_) {
2578            if (fabs(lowerValue-value)<=fabs(upperValue-value)) {
2579              if (upperValue > lowerValue + dualBound_) {
2580                upper_[iSequence]=lowerValue+dualBound_;
2581                setFakeBound(iSequence,ClpSimplexDual::upperFake);
2582                numberFake_++;
2583              }
2584            } else {
2585              if (lowerValue < upperValue - dualBound_) {
2586                lower_[iSequence]=upperValue-dualBound_;
2587                setFakeBound(iSequence,ClpSimplexDual::lowerFake);
2588                numberFake_++;
2589              }
2590            }
2591          } else {
2592            lower_[iSequence]=-0.5*dualBound_;
2593            upper_[iSequence]= 0.5*dualBound_;
2594            setFakeBound(iSequence,ClpSimplexDual::bothFake);
2595            numberFake_++;
2596          }
2597          if (status==atUpperBound)
2598            solution_[iSequence]=upper_[iSequence];
2599          else
2600            solution_[iSequence]=lower_[iSequence];
2601        } else {
2602          // set non basic free variables to fake bounds
2603          // I don't think we should ever get here
2604          CoinAssert(!("should not be here"));
2605          lower_[iSequence]=-0.5*dualBound_;
2606          upper_[iSequence]= 0.5*dualBound_;
2607          setFakeBound(iSequence,ClpSimplexDual::bothFake);
2608          numberFake_++;
2609          setStatus(iSequence,atUpperBound);
2610          solution_[iSequence]=0.5*dualBound_;
2611        }
2612      }
2613    }
2614
2615    return 1;
2616  } else {
2617    // just reset changed ones
2618    if (columnScale_) {
2619      int iSequence;
2620      for (iSequence=0;iSequence<numberColumns_;iSequence++) {
2621        FakeBound fakeStatus = getFakeBound(iSequence);
2622        if (fakeStatus!=noFake) {
2623          if (((int) fakeStatus&1)!=0) {
2624            // lower
2625            double value = columnLower_[iSequence];
2626            if (value>-1.0e30) {
2627              double multiplier = rhsScale_/columnScale_[iSequence];
2628              value *= multiplier;
2629            }
2630            columnLowerWork_[iSequence]=value;
2631          }
2632          if (((int) fakeStatus&2)!=0) {
2633            // upper
2634            double value = columnUpper_[iSequence];
2635            if (value<1.0e30) {
2636              double multiplier = rhsScale_/columnScale_[iSequence];
2637              value *= multiplier;
2638            }
2639            columnUpperWork_[iSequence]=value;
2640          }
2641        }
2642      }
2643      for (iSequence=0;iSequence<numberRows_;iSequence++) {
2644        FakeBound fakeStatus = getFakeBound(iSequence+numberColumns_);
2645        if (fakeStatus!=noFake) {
2646          if (((int) fakeStatus&1)!=0) {
2647            // lower
2648            double value = rowLower_[iSequence];
2649            if (value>-1.0e30) {
2650              double multiplier = rhsScale_*rowScale_[iSequence];
2651              value *= multiplier;
2652            }
2653            rowLowerWork_[iSequence]=value;
2654          }
2655          if (((int) fakeStatus&2)!=0) {
2656            // upper
2657            double value = rowUpper_[iSequence];
2658            if (value<1.0e30) {
2659              double multiplier = rhsScale_*rowScale_[iSequence];
2660              value *= multiplier;
2661            }
2662            rowUpperWork_[iSequence]=value;
2663          }
2664        }
2665      }
2666    } else {
2667      int iSequence;
2668      for (iSequence=0;iSequence<numberColumns_;iSequence++) {
2669        FakeBound fakeStatus = getFakeBound(iSequence);
2670        if (((int) fakeStatus&1)!=0) {
2671          // lower
2672          columnLowerWork_[iSequence]=columnLower_[iSequence];
2673        }
2674        if (((int) fakeStatus&2)!=0) {
2675          // upper
2676          columnUpperWork_[iSequence]=columnUpper_[iSequence];
2677        }
2678      }
2679      for (iSequence=0;iSequence<numberRows_;iSequence++) {
2680        FakeBound fakeStatus = getFakeBound(iSequence+numberColumns_);
2681        if (((int) fakeStatus&1)!=0) {
2682          // lower
2683          rowLowerWork_[iSequence]=rowLower_[iSequence];
2684        }
2685        if (((int) fakeStatus&2)!=0) {
2686          // upper
2687          rowUpperWork_[iSequence]=rowUpper_[iSequence];
2688        }
2689      }
2690    }
2691    return 0;
2692  }
2693}
2694int
2695ClpSimplexDual::dualColumn0(const CoinIndexedVector * rowArray,
2696                           const CoinIndexedVector * columnArray,
2697                           CoinIndexedVector * spareArray,
2698                           double acceptablePivot,
2699                           double & upperReturn, double &bestReturn,double & badFree)
2700{
2701  // do first pass to get possibles
2702  double * spare = spareArray->denseVector();
2703  int * index = spareArray->getIndices();
2704  const double * work;
2705  int number;
2706  const int * which;
2707  const double * reducedCost;
2708  // We can also see if infeasible or pivoting on free
2709  double tentativeTheta = 1.0e25;
2710  double upperTheta = 1.0e31;
2711  double freePivot = acceptablePivot;
2712  double bestPossible=0.0;
2713  int numberRemaining=0;
2714  int i;
2716  for (int iSection=0;iSection<2;iSection++) {
2717
2719
2720    if (!iSection) {
2721      work = rowArray->denseVector();
2722      number = rowArray->getNumElements();
2723      which = rowArray->getIndices();
2724      reducedCost = rowReducedCost_;
2726    } else {
2727      work = columnArray->denseVector();
2728      number = columnArray->getNumElements();
2729      which = columnArray->getIndices();
2730      reducedCost = reducedCostWork_;
2732    }
2733
2734    for (i=0;i<number;i++) {
2735      int iSequence = which[i];
2736      double alpha;
2737      double oldValue;
2738      double value;
2739      bool keep;
2740
2742
2743      case basic:
2744      case ClpSimplex::isFixed:
2745        break;
2746      case isFree:
2747      case superBasic:
2748        alpha = work[i];
2749        bestPossible = CoinMax(bestPossible,fabs(alpha));
2750        oldValue = reducedCost[iSequence];
2751        // If free has to be very large - should come in via dualRow
2753        //break;
2754        if (oldValue>dualTolerance_) {
2755          keep = true;
2756        } else if (oldValue<-dualTolerance_) {
2757          keep = true;
2758        } else {
2759          if (fabs(alpha)>CoinMax(10.0*acceptablePivot,1.0e-5)) {
2760            keep = true;
2761          } else {
2762            keep=false;
2764          }
2765        }
2766        if (keep) {
2767          // free - choose largest
2768          if (fabs(alpha)>freePivot) {
2769            freePivot=fabs(alpha);
2770            sequenceIn_ = iSequence + addSequence;
2771            theta_=oldValue/alpha;
2772            alpha_=alpha;
2773          }
2774        }
2775        break;
2776      case atUpperBound:
2777        alpha = work[i];
2778        oldValue = reducedCost[iSequence];
2779        value = oldValue-tentativeTheta*alpha;
2780        //assert (oldValue<=dualTolerance_*1.0001);
2781        if (value>dualTolerance_) {
2782          bestPossible = CoinMax(bestPossible,-alpha);
2783          value = oldValue-upperTheta*alpha;
2784          if (value>dualTolerance_ && -alpha>=acceptablePivot)
2785            upperTheta = (oldValue-dualTolerance_)/alpha;
2787          spare[numberRemaining]=alpha;
2789        }
2790        break;
2791      case atLowerBound:
2792        alpha = work[i];
2793        oldValue = reducedCost[iSequence];
2794        value = oldValue-tentativeTheta*alpha;
2795        //assert (oldValue>=-dualTolerance_*1.0001);
2796        if (value<-dualTolerance_) {
2797          bestPossible = CoinMax(bestPossible,alpha);
2798          value = oldValue-upperTheta*alpha;
2799          if (value<-dualTolerance_ && alpha>=acceptablePivot)
2800            upperTheta = (oldValue+dualTolerance_)/alpha;
2802          spare[numberRemaining]=alpha;
2804        }
2805        break;
2806      }
2807    }
2808  }
2809  upperReturn = upperTheta;
2810  bestReturn = bestPossible;
2811  return numberRemaining;
2812}
2813/*
2814   Row array has row part of pivot row (as duals so sign may be switched)
2815   Column array has column part.
2816   This chooses pivot column.
2817   Spare array will be needed when we start getting clever.
2818   We will check for basic so spare array will never overflow.
2819   If necessary will modify costs
2820*/
2821double
2822ClpSimplexDual::dualColumn(CoinIndexedVector * rowArray,
2823                           CoinIndexedVector * columnArray,
2824                           CoinIndexedVector * spareArray,
2825                           CoinIndexedVector * spareArray2,
2826                           double acceptablePivot,
2827                           CoinBigIndex * dubiousWeights)
2828{
2829  int numberPossiblySwapped=0;
2830  int numberRemaining=0;
2831
2832  double totalThru=0.0; // for when variables flip
2833  //double saveAcceptable=acceptablePivot;
2834  //acceptablePivot=1.0e-9;
2835
2836  double bestEverPivot=acceptablePivot;
2837  int lastSequence = -1;
2838  double lastPivot=0.0;
2839  double upperTheta;
2840  double newTolerance = dualTolerance_;
2841  //newTolerance = dualTolerance_+1.0e-6*dblParam_[ClpDualTolerance];
2842  // will we need to increase tolerance
2843  bool thisIncrease=false;
2844  // If we think we need to modify costs (not if something from broad sweep)
2845  bool modifyCosts=false;
2846  // Increase in objective due to swapping bounds (may be negative)
2847  double increaseInObjective=0.0;
2848
2849  // use spareArrays to put ones looked at in
2850  // we are going to flip flop between
2851  int iFlip = 0;
2852  // Possible list of pivots
2853  int interesting[2];
2854  // where possible swapped ones are
2855  int swapped[2];
2856  // for zeroing out arrays after
2857  int marker[2][2];
2858  // pivot elements
2859  double * array[2], * spare, * spare2;
2860  // indices
2861  int * indices[2], * index, * index2;
2862  spareArray2->clear();
2863  array[0] = spareArray->denseVector();
2864  indices[0] = spareArray->getIndices();
2865  spare = array[0];
2866  index = indices[0];
2867  array[1] = spareArray2->denseVector();
2868  indices[1] = spareArray2->getIndices();
2869  int i;
2870
2871  // initialize lists
2872  for (i=0;i<2;i++) {
2873    interesting[i]=0;
2874    swapped[i]=numberColumns_;
2875    marker[i][0]=0;
2876    marker[i][1]=numberColumns_;
2877  }
2878  /*
2879    First we get a list of possible pivots.  We can also see if the
2880    problem looks infeasible or whether we want to pivot in free variable.
2881    This may make objective go backwards but can only happen a finite
2882    number of times and I do want free variables basic.
2883
2884    Then we flip back and forth.  At the start of each iteration
2885    interesting[iFlip] should have possible candidates and swapped[iFlip]
2886    will have pivots if we decide to take a previous pivot.
2887    At end of each iteration interesting[1-iFlip] should have
2888    candidates if we go through this theta and swapped[1-iFlip]
2889    pivots if we don't go through.
2890
2891    At first we increase theta and see what happens.  We start
2892    theta at a reasonable guess.  If in right area then we do bit by bit.
2893
2894   */
2895
2896  // do first pass to get possibles
2897  upperTheta = 1.0e31;
2898  double bestPossible=0.0;
2900  if (spareIntArray_[0]!=-1) {
2901    numberRemaining = dualColumn0(rowArray,columnArray,spareArray,
2903  } else {
2905    numberRemaining = spareArray->getNumElements();
2906    spareArray->setNumElements(0);
2907    upperTheta = spareDoubleArray_[0];
2908    bestPossible = spareDoubleArray_[1];
2909    theta_ = spareDoubleArray_[2];
2910    alpha_ = spareDoubleArray_[3];
2911    sequenceIn_ = spareIntArray_[1];
2912  }
2913  // switch off
2914  spareIntArray_[0]=0;
2915  // We can also see if infeasible or pivoting on free
2916  double tentativeTheta = 1.0e25;
2917  interesting[0]=numberRemaining;
2918  marker[0][0] = numberRemaining;
2919
2920  if (!numberRemaining&&sequenceIn_<0)
2921    return 0.0; // Looks infeasible
2922
2923  // If sum of bad small pivots too much
2924#define MORE_CAREFUL
2925#ifdef MORE_CAREFUL
2927#endif
2928  if (sequenceIn_>=0) {
2929    // free variable - always choose
2930  } else {
2931
2932    theta_=1.0e50;
2933    // now flip flop between spare arrays until reasonable theta
2934    tentativeTheta = CoinMax(10.0*upperTheta,1.0e-7);
2935
2936    // loops increasing tentative theta until can't go through
2937
2938    while (tentativeTheta < 1.0e22) {
2939      double thruThis = 0.0;
2940
2941      double bestPivot=acceptablePivot;
2942      int bestSequence=-1;
2943
2944      numberPossiblySwapped = numberColumns_;
2945      numberRemaining = 0;
2946
2947      upperTheta = 1.0e50;
2948
2949      spare = array[iFlip];
2950      index = indices[iFlip];
2951      spare2 = array[1-iFlip];
2952      index2 = indices[1-iFlip];
2953
2954      // try 3 different ways
2955      // 1 bias increase by ones with slightly wrong djs
2956      // 2 bias by all
2957      // 3 bias by all - tolerance
2958#define TRYBIAS 3
2959
2960
2961      double increaseInThis=0.0; //objective increase in this loop
2962
2963      for (i=0;i<interesting[iFlip];i++) {
2964        int iSequence = index[i];
2965        double alpha = spare[i];
2966        double oldValue = dj_[iSequence];
2967        double value = oldValue-tentativeTheta*alpha;
2968
2969        if (alpha < 0.0) {
2970          //at upper bound
2971          if (value>newTolerance) {
2972            double range = upper_[iSequence] - lower_[iSequence];
2973            thruThis -= range*alpha;
2974#if TRYBIAS==1
2975            if (oldValue>0.0)
2976              increaseInThis -= oldValue*range;
2977#elif TRYBIAS==2
2978            increaseInThis -= oldValue*range;
2979#else
2980            increaseInThis -= (oldValue+dualTolerance_)*range;
2981#endif
2982            // goes on swapped list (also means candidates if too many)
2983            spare2[--numberPossiblySwapped]=alpha;
2984            index2[numberPossiblySwapped]=iSequence;
2985            if (fabs(alpha)>bestPivot) {
2986              bestPivot=fabs(alpha);
2987              bestSequence=numberPossiblySwapped;
2988            }
2989          } else {
2990            value = oldValue-upperTheta*alpha;
2991            if (value>newTolerance && -alpha>=acceptablePivot)
2992              upperTheta = (oldValue-newTolerance)/alpha;
2993            spare2[numberRemaining]=alpha;
2994            index2[numberRemaining++]=iSequence;
2995          }
2996        } else {
2997          // at lower bound
2998          if (value<-newTolerance) {
2999            double range = upper_[iSequence] - lower_[iSequence];
3000            thruThis += range*alpha;
3001            //?? is this correct - and should we look at good ones
3002#if TRYBIAS==1
3003            if (oldValue<0.0)
3004              increaseInThis += oldValue*range;
3005#elif TRYBIAS==2
3006            increaseInThis += oldValue*range;
3007#else
3008            increaseInThis += (oldValue-dualTolerance_)*range;
3009#endif
3010            // goes on swapped list (also means candidates if too many)
3011            spare2[--numberPossiblySwapped]=alpha;
3012            index2[numberPossiblySwapped]=iSequence;
3013            if (fabs(alpha)>bestPivot) {
3014              bestPivot=fabs(alpha);
3015              bestSequence=numberPossiblySwapped;
3016            }
3017          } else {
3018            value = oldValue-upperTheta*alpha;
3019            if (value<-newTolerance && alpha>=acceptablePivot)
3020              upperTheta = (oldValue+newTolerance)/alpha;
3021            spare2[numberRemaining]=alpha;
3022            index2[numberRemaining++]=iSequence;
3023          }
3024        }
3025      }
3026      swapped[1-iFlip]=numberPossiblySwapped;
3027      interesting[1-iFlip]=numberRemaining;
3028      marker[1-iFlip][0]= CoinMax(marker[1-iFlip][0],numberRemaining);
3029      marker[1-iFlip][1]= CoinMin(marker[1-iFlip][1],numberPossiblySwapped);
3030
3031      if (totalThru+thruThis>=fabs(dualOut_)||
3032          increaseInObjective+increaseInThis<0.0) {
3033        // We should be pivoting in this batch
3034        // so compress down to this lot
3035        numberRemaining=0;
3036        for (i=numberColumns_-1;i>=swapped[1-iFlip];i--) {
3037          spare[numberRemaining]=spare2[i];
3038          index[numberRemaining++]=index2[i];
3039        }
3040        interesting[iFlip]=numberRemaining;
3041        int iTry;
3042#define MAXTRY 100
3043        // first get ratio with tolerance
3044        for (iTry=0;iTry<MAXTRY;iTry++) {
3045
3046          upperTheta=1.0e50;
3047          numberPossiblySwapped = numberColumns_;
3048          numberRemaining = 0;
3049
3050          increaseInThis=0.0; //objective increase in this loop
3051
3052          thruThis=0.0;
3053
3054          spare = array[iFlip];
3055          index = indices[iFlip];
3056          spare2 = array[1-iFlip];
3057          index2 = indices[1-iFlip];
3058          for (i=0;i<interesting[iFlip];i++) {
3059            int iSequence=index[i];
3060            double alpha=spare[i];
3061            double oldValue = dj_[iSequence];
3062            double value = oldValue-upperTheta*alpha;
3063
3064            if (alpha < 0.0) {
3065              //at upper bound
3066              if (value>newTolerance) {
3067                if (-alpha>=acceptablePivot) {
3068                  upperTheta = (oldValue-newTolerance)/alpha;
3069                }
3070              }
3071            } else {
3072              // at lower bound
3073              if (value<-newTolerance) {
3074                if (alpha>=acceptablePivot) {
3075                  upperTheta = (oldValue+newTolerance)/alpha;
3076                }
3077              }
3078            }
3079          }
3080          bestPivot=acceptablePivot;
3081          sequenceIn_=-1;
3082#ifdef DUBIOUS_WEIGHTS
3083          double bestWeight=COIN_DBL_MAX;
3084#endif
3085          double largestPivot=acceptablePivot;
3086          // now choose largest and sum all ones which will go through
3087          //printf("XX it %d number %d\n",numberIterations_,interesting[iFlip]);
3088  // Sum of bad small pivots
3089#ifdef MORE_CAREFUL
3092#endif
3093          // Make sure upperTheta will work (-O2 and above gives problems)
3094          upperTheta *= 1.0000000001;
3095          for (i=0;i<interesting[iFlip];i++) {
3096            int iSequence=index[i];
3097            double alpha=spare[i];
3098            double value = dj_[iSequence]-upperTheta*alpha;
3100
3102
3103            if (alpha < 0.0) {
3104              //at upper bound
3105              if (value>=0.0) {
3107#if TRYBIAS==1
3109#elif TRYBIAS==2
3111#else
3113#endif
3114              }
3115            } else {
3116              // at lower bound
3117              if (value<=0.0) {
3119#if TRYBIAS==1
3121#elif TRYBIAS==2
3123#else
3125#endif
3126              }
3127            }
3129              // add to list of remaining
3130              spare2[numberRemaining]=alpha;
3131              index2[numberRemaining++]=iSequence;
3132            } else {
3133              // add to list of swapped
3134              spare2[--numberPossiblySwapped]=alpha;
3135              index2[numberPossiblySwapped]=iSequence;
3136              // select if largest pivot
3137              bool take=false;
3138              double absAlpha = fabs(alpha);
3139#ifdef DUBIOUS_WEIGHTS
3140              // User could do anything to break ties here
3141              double weight;
3142              if (dubiousWeights)
3143                weight=dubiousWeights[iSequence];
3144              else
3145                weight=1.0;
3146              weight += randomNumberGenerator_.randomDouble()*1.0e-2;
3147              if (absAlpha>2.0*bestPivot) {
3148                take=true;
3149              } else if (absAlpha>largestPivot) {
3150                // could multiply absAlpha and weight
3151                if (weight*bestPivot<bestWeight*absAlpha)
3152                  take=true;
3153              }
3154#else
3155              if (absAlpha>bestPivot)
3156                take=true;
3157#endif
3158#ifdef MORE_CAREFUL
3159              if (absAlpha<acceptablePivot&&upperTheta<1.0e20) {
3160                if (alpha < 0.0) {
3161                  //at upper bound
3162                  if (value>dualTolerance_) {
3163                    double gap=upper_[iSequence]-lower_[iSequence];
3164                    if (gap<1.0e20)
3166                    else
3168                    //printf("bad %d alpha %g dj at upper %g\n",
3169                    //     iSequence,alpha,value);
3170                  }
3171                } else {
3172                  //at lower bound
3173                  if (value<-dualTolerance_) {
3174                    double gap=upper_[iSequence]-lower_[iSequence];
3175                    if (gap<1.0e20)
3177                    else
3179                    //printf("bad %d alpha %g dj at lower %g\n",
3180                    //     iSequence,alpha,value);
3181                  }
3182                }
3183              }
3184#endif
3185#ifdef FORCE_FOLLOW
3186              if (iSequence==force_in) {
3187                printf("taking %d - alpha %g best %g\n",force_in,absAlpha,largestPivot);
3188                take=true;
3189              }
3190#endif
3191              if (take) {
3192                sequenceIn_ = numberPossiblySwapped;
3193                bestPivot =  absAlpha;
3194                theta_ = dj_[iSequence]/alpha;
3195                largestPivot = CoinMax(largestPivot,0.5*bestPivot);
3196#ifdef DUBIOUS_WEIGHTS
3197                bestWeight = weight;
3198#endif
3199                //printf(" taken seq %d alpha %g weight %d\n",
3200                //   iSequence,absAlpha,dubiousWeights[iSequence]);
3201              } else {
3202                //printf(" not taken seq %d alpha %g weight %d\n",
3203                //   iSequence,absAlpha,dubiousWeights[iSequence]);
3204              }
3205              double range = upper_[iSequence] - lower_[iSequence];
3206              thruThis += range*fabs(alpha);
3208            }
3209          }
3210#ifdef MORE_CAREFUL
3211          // If we have done pivots and things look bad set alpha_ 0.0 to force factorization
3213            if (handler_->logLevel()>1)
3214            printf("maybe forcing re-factorization - sum %g  %d pivots\n",sumBadPivots,
3215                   factorization_->pivots());
3217          }
3218#endif
3219          swapped[1-iFlip]=numberPossiblySwapped;
3220          interesting[1-iFlip]=numberRemaining;
3221          marker[1-iFlip][0]= CoinMax(marker[1-iFlip][0],numberRemaining);
3222          marker[1-iFlip][1]= CoinMin(marker[1-iFlip][1],numberPossiblySwapped);
3223          // If we stop now this will be increase in objective (I think)
3224          double increase = (fabs(dualOut_)-totalThru)*theta_;
3225          increase += increaseInObjective;
3226          if (theta_<0.0)
3227            thruThis += fabs(dualOut_); // force using this one
3228          if (increaseInObjective<0.0&&increase<0.0&&lastSequence>=0) {
3229            // back
3230            // We may need to be more careful - we could do by
3231            // switch so we always do fine grained?
3232            bestPivot=0.0;
3233          } else {
3235            totalThru += thruThis;
3236            increaseInObjective += increaseInThis;
3237          }
3238          if (bestPivot<0.1*bestEverPivot&&
3239              bestEverPivot>1.0e-6&&
3240              (bestPivot<1.0e-3||totalThru*2.0>fabs(dualOut_))) {
3241            // back to previous one
3242            sequenceIn_=lastSequence;
3243            // swap regions
3244            iFlip = 1-iFlip;
3245            break;
3246          } else if (sequenceIn_==-1&&upperTheta>largeValue_) {
3247            if (lastPivot>acceptablePivot) {
3248              // back to previous one
3249              sequenceIn_=lastSequence;
3250              // swap regions
3251              iFlip = 1-iFlip;
3252            } else {
3253              // can only get here if all pivots too small
3254            }
3255            break;
3256          } else if (totalThru>=fabs(dualOut_)) {
3257            modifyCosts=true; // fine grain - we can modify costs
3258            break; // no point trying another loop
3259          } else {
3260            lastSequence=sequenceIn_;
3261            if (bestPivot>bestEverPivot)
3262              bestEverPivot=bestPivot;
3263            iFlip = 1 -iFlip;
3264            modifyCosts=true; // fine grain - we can modify costs
3265          }
3266        }
3267        if (iTry==MAXTRY)
3268          iFlip = 1-iFlip; // flip back
3269        break;
3270      } else {
3271        // skip this lot
3272        if (bestPivot>1.0e-3||bestPivot>bestEverPivot) {
3273          bestEverPivot=bestPivot;
3274          lastSequence=bestSequence;
3275        } else {
3276          // keep old swapped
3277          CoinMemcpyN(array[iFlip]+swapped[iFlip],
3278                 numberColumns_-swapped[iFlip],array[1-iFlip]+swapped[iFlip]);
3279          CoinMemcpyN(indices[iFlip]+swapped[iFlip],
3280                 numberColumns_-swapped[iFlip],indices[1-iFlip]+swapped[iFlip]);
3281          marker[1-iFlip][1] = CoinMin(marker[1-iFlip][1],swapped[iFlip]);
3282          swapped[1-iFlip]=swapped[iFlip];
3283        }
3284        increaseInObjective += increaseInThis;
3285        iFlip = 1 - iFlip; // swap regions
3286        tentativeTheta = 2.0*upperTheta;
3287        totalThru += thruThis;
3288      }
3289    }
3290
3291    // can get here without sequenceIn_ set but with lastSequence
3292    if (sequenceIn_<0&&lastSequence>=0) {
3293      // back to previous one
3294      sequenceIn_=lastSequence;
3295      // swap regions
3296      iFlip = 1-iFlip;
3297    }
3298
3299#define MINIMUMTHETA 1.0e-18
3300    // Movement should be minimum for anti-degeneracy - unless
3301    // fixed variable out
3302    double minimumTheta;
3303    if (upperOut_>lowerOut_)
3304      minimumTheta=MINIMUMTHETA;
3305    else
3306      minimumTheta=0.0;
3307    if (sequenceIn_>=0) {
3308      // at this stage sequenceIn_ is just pointer into index array
3309      // flip just so we can use iFlip
3310      iFlip = 1 -iFlip;
3311      spare = array[iFlip];
3312      index = indices[iFlip];
3313      double oldValue;
3314      alpha_ = spare[sequenceIn_];
3315      sequenceIn_ = indices[iFlip][sequenceIn_];
3316      oldValue = dj_[sequenceIn_];
3317      theta_ = CoinMax(oldValue/alpha_,0.0);
3318      if (theta_<minimumTheta&&fabs(alpha_)<1.0e5&&1) {
3319        // can't pivot to zero
3320#if 0
3321        if (oldValue-minimumTheta*alpha_>=-dualTolerance_) {
3322          theta_=minimumTheta;
3323        } else if (oldValue-minimumTheta*alpha_>=-newTolerance) {
3324          theta_=minimumTheta;
3325          thisIncrease=true;
3326        } else {
3327          theta_=CoinMax((oldValue+newTolerance)/alpha_,0.0);
3328          thisIncrease=true;
3329        }
3330#else
3331        theta_=minimumTheta;
3332#endif
3333      }
3334      // may need to adjust costs so all dual feasible AND pivoted is exactly 0
3335      //int costOffset = numberRows_+numberColumns_;
3336      if (modifyCosts) {
3337        int i;
3338        for (i=numberColumns_-1;i>=swapped[iFlip];i--) {
3339          int iSequence=index[i];
3340          double alpha=spare[i];
3341          double value = dj_[iSequence]-theta_*alpha;
3342
3343          // can't be free here
3344
3345          if (alpha < 0.0) {
3346            //at upper bound
3347            if (value>dualTolerance_) {
3348              thisIncrease=true;
3349#define MODIFYCOST 2
3350#if MODIFYCOST
3351              // modify cost to hit new tolerance
3352              double modification = alpha*theta_-dj_[iSequence]
3353                +newTolerance;
3354              if ((specialOptions_&(2048+4096+16384))!=0) {
3355                if ((specialOptions_&16384)!=0) {
3356                  if (fabs(modification)<1.0e-8)
3357                    modification=0.0;
3358                } else if ((specialOptions_&2048)!=0) {
3359                  if (fabs(modification)<1.0e-10)
3360                    modification=0.0;
3361                } else {
3362                  if (fabs(modification)<1.0e-12)
3363                    modification=0.0;
3364                }
3365              }
3366              dj_[iSequence] += modification;
3367              cost_[iSequence] +=  modification;
3368              if (modification)
3369                numberChanged_ ++; // Say changed costs
3370              //cost_[iSequence+costOffset] += modification; // save change
3371#endif
3372            }
3373          } else {
3374            // at lower bound
3375            if (-value>dualTolerance_) {
3376              thisIncrease=true;
3377#if MODIFYCOST
3378              // modify cost to hit new tolerance
3379              double modification = alpha*theta_-dj_[iSequence]
3380                -newTolerance;
3381              //modification = CoinMax(modification,-dualTolerance_);
3382              //assert (fabs(modification)<1.0e-7);
3383              if ((specialOptions_&(2048+4096))!=0) {
3384                if ((specialOptions_&2048)!=0) {
3385                  if (fabs(modification)<1.0e-10)
3386                    modification=0.0;
3387                } else {
3388                  if (fabs(modification)<1.0e-12)
3389                    modification=0.0;
3390                }
3391              }
3392              dj_[iSequence] += modification;
3393              cost_[iSequence] +=  modification;
3394              if (modification)
3395                numberChanged_ ++; // Say changed costs
3396              //cost_[iSequence+costOffset] += modification; // save change
3397#endif
3398            }
3399          }
3400        }
3401      }
3402    }
3403  }
3404
3405  if (sequenceIn_>=0) {
3406#ifdef MORE_CAREFUL
3407    // If we have done pivots and things look bad set alpha_ 0.0 to force factorization
3409      if (handler_->logLevel()>1)
3410        printf("forcing re-factorization\n");
3411      alpha_=0.0;
3412    }
3414      if (handler_->logLevel()>1)
3415        printf("forcing re-factorizationon free\n");
3416      alpha_=0.0;
3417    }
3418#endif
3419    lowerIn_ = lower_[sequenceIn_];
3420    upperIn_ = upper_[sequenceIn_];
3421    valueIn_ = solution_[sequenceIn_];
3422    dualIn_ = dj_[sequenceIn_];
3423
3424    if (numberTimesOptimal_) {
3425      // can we adjust cost back closer to original
3427    }
3428#if MODIFYCOST>1
3429    // modify cost to hit zero exactly
3430    // so (dualIn_+modification)==theta_*alpha_
3431    double modification = theta_*alpha_-dualIn_;
3432    // But should not move objectivetoo much ??
3433#define DONT_MOVE_OBJECTIVE
3434#ifdef DONT_MOVE_OBJECTIVE
3435    double moveObjective = fabs(modification*solution_[sequenceIn_]);
3436    double smallMove = CoinMax(fabs(objectiveValue_),1.0e-3);
3437    if (moveObjective>smallMove) {
3438      if (handler_->logLevel()>1)
3439        printf("would move objective by %g - original mod %g sol value %g\n",moveObjective,
3440               modification,solution_[sequenceIn_]);
3441      modification *= smallMove/moveObjective;
3442    }
3443#endif
3444    if ((specialOptions_&(2048+4096))!=0) {
3445      if ((specialOptions_&16384)!=0) {
3446        // in fast dual
3447        if (fabs(modification)<1.0e-7)
3448          modification=0.0;
3449      } else if ((specialOptions_&2048)!=0) {
3450        if (fabs(modification)<1.0e-10)
3451          modification=0.0;
3452      } else {
3453        if (fabs(modification)<1.0e-12)
3454          modification=0.0;
3455      }
3456    }
3457    dualIn_ += modification;
3458    dj_[sequenceIn_]=dualIn_;
3459    cost_[sequenceIn_] += modification;
3460    if (modification)
3461      numberChanged_ ++; // Say changed costs
3462    //int costOffset = numberRows_+numberColumns_;
3463    //cost_[sequenceIn_+costOffset] += modification; // save change
3464    //assert (fabs(modification)<1.0e-6);
3465#ifdef CLP_DEBUG
3466    if ((handler_->logLevel()&32)&&fabs(modification)>1.0e-15)
3467      printf("exact %d new cost %g, change %g\n",sequenceIn_,
3468             cost_[sequenceIn_],modification);
3469#endif
3470#endif
3471
3472    if (alpha_<0.0) {
3473      // as if from upper bound
3474      directionIn_=-1;
3475      upperIn_=valueIn_;
3476    } else {
3477      // as if from lower bound
3478      directionIn_=1;
3479      lowerIn_=valueIn_;
3480    }
3481  }
3482  //if (thisIncrease)
3483  //dualTolerance_+= 1.0e-6*dblParam_[ClpDualTolerance];
3484
3485  // clear arrays
3486
3487  for (i=0;i<2;i++) {
3488    CoinZeroN(array[i],marker[i][0]);
3489    CoinZeroN(array[i]+marker[i][1],numberColumns_-marker[i][1]);
3490  }
3491  return bestPossible;
3492}
3493#ifdef CLP_ALL_ONE_FILE
3494#undef MAXTRY
3495#endif
3496/* Checks if tentative optimal actually means unbounded
3497   Returns -3 if not, 2 if is unbounded */
3498int
3499ClpSimplexDual::checkUnbounded(CoinIndexedVector * ray,
3500                                   CoinIndexedVector * spare,
3501                                   double changeCost)
3502{
3503  int status=2; // say unbounded
3504  factorization_->updateColumn(spare,ray);
3505  // get reduced cost
3506  int i;
3507  int number=ray->getNumElements();
3508  int * index = ray->getIndices();
3509  double * array = ray->denseVector();
3510  for (i=0;i<number;i++) {
3511    int iRow=index[i];
3512    int iPivot=pivotVariable_[iRow];
3513    changeCost -= cost(iPivot)*array[iRow];
3514  }
3515  double way;
3516  if (changeCost>0.0) {
3517    //try going down
3518    way=1.0;
3519  } else if (changeCost<0.0) {
3520    //try going up
3521    way=-1.0;
3522  } else {
3523#ifdef CLP_DEBUG
3524    printf("can't decide on up or down\n");
3525#endif
3526    way=0.0;
3527    status=-3;
3528  }
3529  double movement=1.0e10*way; // some largish number
3530  double zeroTolerance = 1.0e-14*dualBound_;
3531  for (i=0;i<number;i++) {
3532    int iRow=index[i];
3533    int iPivot=pivotVariable_[iRow];
3534    double arrayValue = array[iRow];
3535    if (fabs(arrayValue)<zeroTolerance)
3536      arrayValue=0.0;
3537    double newValue=solution(iPivot)+movement*arrayValue;
3538    if (newValue>upper(iPivot)+primalTolerance_||
3539        newValue<lower(iPivot)-primalTolerance_)
3540      status=-3; // not unbounded
3541  }
3542  if (status==2) {
3543    // create ray
3544    delete [] ray_;
3545    ray_ = new double [numberColumns_];
3546    CoinZeroN(ray_,numberColumns_);
3547    for (i=0;i<number;i++) {
3548      int iRow=index[i];
3549      int iPivot=pivotVariable_[iRow];
3550      double arrayValue = array[iRow];
3551      if (iPivot<numberColumns_&&fabs(arrayValue)>=zeroTolerance)
3552        ray_[iPivot] = way* array[iRow];
3553    }
3554  }
3555  ray->clear();
3556  return status;
3557}
3558//static int count_alpha=0;
3559/* Checks if finished.  Updates status */
3560void
3561ClpSimplexDual::statusOfProblemInDual(int & lastCleaned,int type,
3562                                      double * givenDuals, ClpDataSave & saveData,
3563                                      int ifValuesPass)
3564{
3565#ifdef CLP_INVESTIGATE
3566  if (z_thinks>0&&z_thinks<2)
3567    z_thinks+=2;
3568#endif
3569  // If lots of iterations then adjust costs if large ones
3570  if (numberIterations_>4*(numberRows_+numberColumns_)&&objectiveScale_==1.0) {
3571    double largest=0.0;
3572    for (int i=0;i<numberRows_;i++) {
3573      int iColumn = pivotVariable_[i];
3574      largest = CoinMax(largest,fabs(cost_[iColumn]));
3575    }
3576    if (largest>1.0e6) {
3577      objectiveScale_ = 1.0e6/largest;
3578      for (int i=0;i<numberRows_+numberColumns_;i++)
3579        cost_[i] *= objectiveScale_;
3580    }
3581  }
3582  int numberPivots = factorization_->pivots();
3583  double realDualInfeasibilities=0.0;
3584  if (type==2) {
3585    if (alphaAccuracy_!=-1.0)
3586      alphaAccuracy_=-2.0;
3587    // trouble - restore solution
3588    CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
3589    CoinMemcpyN(savedSolution_+numberColumns_ ,
3590           numberRows_,rowActivityWork_);
3591    CoinMemcpyN(savedSolution_ ,
3592           numberColumns_,columnActivityWork_);
3593    // restore extra stuff
3594    int dummy;
3595    matrix_->generalExpanded(this,6,dummy);
3596    forceFactorization_=1; // a bit drastic but ..
3597    changeMade_++; // say something changed
3598    // get correct bounds on all variables
3599    resetFakeBounds();
3600  }
3601  int tentativeStatus = problemStatus_;
3602  double changeCost;
3603  bool unflagVariables = true;
3604  bool weightsSaved=false;
3605  int dontFactorizePivots = dontFactorizePivots_;
3606  if (type==3) {
3607    type=1;
3608    dontFactorizePivots=1;
3609  }
3610  if (alphaAccuracy_<0.0||!numberPivots||alphaAccuracy_>1.0e4||numberPivots>20) {
3611    if (problemStatus_>-3||numberPivots>dontFactorizePivots) {
3612      // factorize
3613      // later on we will need to recover from singularities
3614      // also we could skip if first time
3615      // save dual weights
3616      dualRowPivot_->saveWeights(this,1);
3617      weightsSaved=true;
3618      if (type) {
3619        // is factorization okay?
3620        if (internalFactorize(1)) {
3621          // no - restore previous basis
3622          unflagVariables = false;
3623          assert (type==1);
3624          changeMade_++; // say something changed
3625          // Keep any flagged variables
3626          int i;
3627          for (i=0;i<numberRows_+numberColumns_;i++) {
3628            if (flagged(i))
3629              saveStatus_[i] |= 64; //say flagged
3630          }
3631          CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
3632          CoinMemcpyN(savedSolution_+numberColumns_ ,
3633                      numberRows_,rowActivityWork_);
3634          CoinMemcpyN(savedSolution_ ,
3635                      numberColumns_,columnActivityWork_);
3636          // restore extra stuff
3637          int dummy;
3638          matrix_->generalExpanded(this,6,dummy);
3639          // get correct bounds on all variables
3640          resetFakeBounds();
3641          // need to reject something
3642          char x = isColumn(sequenceOut_) ? 'C' :'R';
3643          handler_->message(CLP_SIMPLEX_FLAG,messages_)
3644            <<x<<sequenceWithin(sequenceOut_)
3645            <<CoinMessageEol;
3646#ifdef COIN_DEVELOP
3647          printf("flag d\n");
3648#endif
3649          setFlagged(sequenceOut_);
3651
3652          // Go to safe
3653          factorization_->pivotTolerance(0.99);
3654          forceFactorization_=1; // a bit drastic but ..
3655          type = 2;
3656          //assert (internalFactorize(1)==0);
3657          if (internalFactorize(1)) {
3658            CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
3659            CoinMemcpyN(savedSolution_+numberColumns_ ,
3660                        numberRows_,rowActivityWork_);
3661            CoinMemcpyN(savedSolution_ ,
3662                        numberColumns_,columnActivityWork_);
3663            // restore extra stuff
3664            int dummy;
3665            matrix_->generalExpanded(this,6,dummy);
3666            // debug
3667            int returnCode = internalFactorize(1);
3668            while (returnCode) {
3669              // ouch
3670              // switch off dense
3671              int saveDense = factorization_->denseThreshold();
3672              factorization_->setDenseThreshold(0);
3673              // Go to safe
3674              factorization_->pivotTolerance(0.99);
3675              // make sure will do safe factorization
3676              pivotVariable_[0]=-1;
3677              returnCode=internalFactorize(2);
3678              factorization_->setDenseThreshold(saveDense);
3679            }
3680            // get correct bounds on all variables
3681            resetFakeBounds();
3682          }
3683        }
3684      }
3685      if (problemStatus_!=-4||numberPivots>10)
3686        problemStatus_=-3;
3687    }
3688  } else {
3689    //printf("testing with accuracy of %g and status of %d\n",alphaAccuracy_,problemStatus_);
3690    //count_alpha++;
3691    //if ((count_alpha%5000)==0)
3692    //printf("count alpha %d\n",count_alpha);
3693  }
3694  // at this stage status is -3 or -4 if looks infeasible
3695  // get primal and dual solutions
3696#if 0
3697  {
3698    int numberTotal = numberRows_+numberColumns_;
3699    double * saveSol = CoinCopyOfArray(solution_,numberTotal);
3700    double * saveDj = CoinCopyOfArray(dj_,numberTotal);
3701    double tolerance = type ? 1.0e-4 : 1.0e-8;
3702    // always if values pass
3703    double saveObj=objectiveValue_;
3704    double sumPrimal = sumPrimalInfeasibilities_;
3705    int numberPrimal = numberPrimalInfeasibilities_;
3706    double sumDual = sumDualInfeasibilities_;
3707    int numberDual = numberDualInfeasibilities_;
3708    gutsOfSolution(givenDuals,NULL);
3709    int j;
3710    double largestPrimal = tolerance;
3711    int iPrimal=-1;
3712    for (j=0;j<numberTotal;j++) {
3713      double difference = solution_[j]-saveSol[j];
3714      if (fabs(difference)>largestPrimal) {
3715        iPrimal=j;
3716        largestPrimal=fabs(difference);
3717      }
3718    }
3719    double largestDual = tolerance;
3720    int iDual=-1;
3721    for (j=0;j<numberTotal;j++) {
3722      double difference = dj_[j]-saveDj[j];
3723      if (fabs(difference)>largestDual&&upper_[j]>lower_[j]) {
3724        iDual=j;
3725        largestDual=fabs(difference);
3726      }
3727    }
3728    if (!type) {
3729      if (fabs(saveObj-objectiveValue_)>1.0e-5||
3730          numberPrimal!=numberPrimalInfeasibilities_||numberPrimal!=1||
3731          fabs(sumPrimal-sumPrimalInfeasibilities_)>1.0e-5||iPrimal>=0||
3732          numberDual!=numberDualInfeasibilities_||numberDual!=0||
3733          fabs(sumDual-sumDualInfeasibilities_)>1.0e-5||iDual>=0)
3734        printf("type %d its %d pivots %d primal n(%d,%d) s(%g,%g) diff(%g,%d) dual n(%d,%d) s(%g,%g) diff(%g,%d) obj(%g,%g)\n",
3735               type,numberIterations_,numberPivots,
3736               numberPrimal,numberPrimalInfeasibilities_,sumPrimal,sumPrimalInfeasibilities_,
3737               largestPrimal,iPrimal,
3738               numberDual,numberDualInfeasibilities_,sumDual,sumDualInfeasibilities_,
3739               largestDual,iDual,
3740               saveObj,objectiveValue_);
3741    } else {
3742      if (fabs(saveObj-objectiveValue_)>1.0e-5||
3743          numberPrimalInfeasibilities_||iPrimal>=0||
3744          numberDualInfeasibilities_||iDual>=0)
3745        printf("type %d its %d pivots %d primal n(%d,%d) s(%g,%g) diff(%g,%d) dual n(%d,%d) s(%g,%g) diff(%g,%d) obj(%g,%g)\n",
3746               type,numberIterations_,numberPivots,
3747               numberPrimal,numberPrimalInfeasibilities_,sumPrimal,sumPrimalInfeasibilities_,
3748               largestPrimal,iPrimal,
3749               numberDual,numberDualInfeasibilities_,sumDual,sumDualInfeasibilities_,
3750               largestDual,iDual,
3751               saveObj,objectiveValue_);
3752    }
3753    delete [] saveSol;
3754    delete [] saveDj;
3755  }
3756#else
3757  if (type||ifValuesPass)
3758    gutsOfSolution(givenDuals,NULL);
3759#endif
3760  // If bad accuracy treat as singular
3761  if ((largestPrimalError_>1.0e15||largestDualError_>1.0e15)&&numberIterations_) {
3762    // restore previous basis
3763    unflagVariables = false;
3764    changeMade_++; // say something changed
3765    // Keep any flagged variables
3766    int i;
3767    for (i=0;i<numberRows_+numberColumns_;i++) {
3768      if (flagged(i))
3769        saveStatus_[i] |= 64; //say flagged
3770    }
3771    CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
3772    CoinMemcpyN(savedSolution_+numberColumns_ ,
3773           numberRows_,rowActivityWork_);
3774    CoinMemcpyN(savedSolution_ ,
3775           numberColumns_,columnActivityWork_);
3776    // restore extra stuff
3777    int dummy;
3778    matrix_->generalExpanded(this,6,dummy);
3779    // get correct bounds on all variables
3780    resetFakeBounds();
3781    // need to reject something
3782    char x = isColumn(sequenceOut_) ? 'C' :'R';
3783    handler_->message(CLP_SIMPLEX_FLAG,messages_)
3784      <<x<<sequenceWithin(sequenceOut_)
3785      <<CoinMessageEol;
3786#ifdef COIN_DEVELOP
3787    printf("flag e\n");
3788#endif
3789    setFlagged(sequenceOut_);
3791
3792    // Go to safer
3793    double newTolerance = CoinMin(1.1*factorization_->pivotTolerance(),0.99);
3794    factorization_->pivotTolerance(newTolerance);
3795    forceFactorization_=1; // a bit drastic but ..
3796    if (alphaAccuracy_!=-1.0)
3797      alphaAccuracy_=-2.0;
3798    type = 2;
3799    //assert (internalFactorize(1)==0);
3800    if (internalFactorize(1)) {
3801      CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
3802      CoinMemcpyN(savedSolution_+numberColumns_ ,
3803              numberRows_,rowActivityWork_);
3804      CoinMemcpyN(savedSolution_ ,
3805           numberColumns_,columnActivityWork_);
3806      // restore extra stuff
3807      int dummy;
3808      matrix_->generalExpanded(this,6,dummy);
3809      // debug
3810      int returnCode = internalFactorize(1);
3811      while (returnCode) {
3812        // ouch
3813        // switch off dense
3814        int saveDense = factorization_->denseThreshold();
3815        factorization_->setDenseThreshold(0);
3816        // Go to safe
3817        factorization_->pivotTolerance(0.99);
3818        // make sure will do safe factorization
3819        pivotVariable_[0]=-1;
3820        returnCode=internalFactorize(2);
3821        factorization_->setDenseThreshold(saveDense);
3822      }
3823      // get correct bounds on all variables
3824      resetFakeBounds();
3825    }
3826    // get primal and dual solutions
3827    gutsOfSolution(givenDuals,NULL);
3828  } else if (goodAccuracy()) {
3829    // Can reduce tolerance
3830    double newTolerance = CoinMax(0.99*factorization_->pivotTolerance(),saveData.pivotTolerance_);
3831    factorization_->pivotTolerance(newTolerance);
3832  }
3833  // Double check infeasibility if no action
3834  if (progress_.lastIterationNumber(0)==numberIterations_) {
3835    if (dualRowPivot_->looksOptimal()) {
3836      numberPrimalInfeasibilities_ = 0;
3837      sumPrimalInfeasibilities_ = 0.0;
3838    }
3839#if 1
3840  } else {
3841    double thisObj = objectiveValue_;
3842    double lastObj = progress_.lastObjective(0);
3843    if(!ifValuesPass&&firstFree_<0) {
3844      if (lastObj>thisObj+1.0e-3*CoinMax(fabs(thisObj),fabs(lastObj))+1.0) {
3845        int maxFactor = factorization_->maximumPivots();
3846        if (maxFactor>10&&numberPivots>1) {
3847          //printf("lastobj %g thisobj %g\n",lastObj,thisObj);
3848          //if (forceFactorization_<0)
3849          //forceFactorization_= maxFactor;
3850          //forceFactorization_ = CoinMax(1,(forceFactorization_>>1));
3851          forceFactorization_=1;
3852          //printf("Reducing factorization frequency - bad backwards\n");
3853          unflagVariables = false;
3854          changeMade_++; // say something changed
3855          CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
3856          CoinMemcpyN(savedSolution_+numberColumns_ ,
3857                      numberRows_,rowActivityWork_);
3858          CoinMemcpyN(savedSolution_ ,
3859                      numberColumns_,columnActivityWork_);
3860          // restore extra stuff
3861          int dummy;
3862          matrix_->generalExpanded(this,6,dummy);
3863          // get correct bounds on all variables
3864          resetFakeBounds();
3865          if(factorization_->pivotTolerance()<0.2)
3866            factorization_->pivotTolerance(0.2);
3867          if (alphaAccuracy_!=-1.0)
3868            alphaAccuracy_=-2.0;
3869          if (internalFactorize(1)) {
3870            CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
3871            CoinMemcpyN(savedSolution_+numberColumns_ ,
3872                        numberRows_,rowActivityWork_);
3873            CoinMemcpyN(savedSolution_ ,
3874                        numberColumns_,columnActivityWork_);
3875            // restore extra stuff
3876            int dummy;
3877            matrix_->generalExpanded(this,6,dummy);
3878            // debug
3879            int returnCode = internalFactorize(1);
3880            while (returnCode) {
3881              // ouch
3882              // switch off dense
3883              int saveDense = factorization_->denseThreshold();
3884              factorization_->setDenseThreshold(0);
3885              // Go to safe
3886              factorization_->pivotTolerance(0.99);
3887              // make sure will do safe factorization
3888              pivotVariable_[0]=-1;
3889              returnCode=internalFactorize(2);
3890              factorization_->setDenseThreshold(saveDense);
3891            }
3892            resetFakeBounds();
3893          }
3894          type = 2; // so will restore weights
3895          // get primal and dual solutions
3896          gutsOfSolution(givenDuals,NULL);
3897        }
3898      } else if (lastObj<thisObj-1.0e-5*CoinMax(fabs(thisObj),fabs(lastObj))-1.0e-3) {
3899        numberTimesOptimal_=0;
3900      }
3901    }
3902#endif
3903  }
3904  // Up tolerance if looks a bit odd
3905  if (numberIterations_>CoinMax(1000,numberRows_>>4)&&(specialOptions_&64)!=0) {
3906    if (sumPrimalInfeasibilities_&&sumPrimalInfeasibilities_<1.0e5) {
3907      int backIteration = progress_.lastIterationNumber(CLP_PROGRESS-1);
3908      if (backIteration>0&&numberIterations_-backIteration<9*CLP_PROGRESS) {
3909        if (factorization_->pivotTolerance()<0.9) {
3910          // up tolerance
3911          factorization_->pivotTolerance(CoinMin(factorization_->pivotTolerance()*1.05+0.02,0.91));
3912          //printf("tol now %g\n",factorization_->pivotTolerance());
3913          progress_.clearIterationNumbers();
3914        }
3915      }
3916    }
3917  }
3918  // Check if looping
3919  int loop;
3920  if (!givenDuals&&type!=2)
3921    loop = progress_.looping();
3922  else
3923    loop=-1;
3924  int situationChanged=0;
3925  if (loop>=0) {
3926    problemStatus_ = loop; //exit if in loop
3927    if (!problemStatus_) {
3928      // declaring victory
3929      numberPrimalInfeasibilities_ = 0;
3930      sumPrimalInfeasibilities_ = 0.0;
3931    } else {
3932      problemStatus_ = 10; // instead - try other algorithm
3933#if COIN_DEVELOP>2
3934      printf("returning at %d\n",__LINE__);
3935#endif
3936    }
3937    return;
3938  } else if (loop<-1) {
3939    // something may have changed
3940    gutsOfSolution(NULL,NULL);
3941    situationChanged=1;
3942  }
3944  if((progressFlag_&2)!=0) {
3945    situationChanged=2;
3946  }
3947  progressFlag_ = 0; //reset progress flag
3948  if (handler_->detail(CLP_SIMPLEX_STATUS,messages_)<100) {
3949    handler_->message(CLP_SIMPLEX_STATUS,messages_)
3950      <<numberIterations_<<objectiveValue();
3951    handler_->printing(sumPrimalInfeasibilities_>0.0)
3952      <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
3953    handler_->printing(sumDualInfeasibilities_>0.0)
3954      <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
3955    handler_->printing(numberDualInfeasibilitiesWithoutFree_
3956                       <numberDualInfeasibilities_)
3957                         <<numberDualInfeasibilitiesWithoutFree_;
3958    handler_->message()<<CoinMessageEol;
3959  }
3960  realDualInfeasibilities=sumDualInfeasibilities_;
3961  double saveTolerance =dualTolerance_;
3962  // If we need to carry on cleaning variables
3963  if (!numberPrimalInfeasibilities_&&(specialOptions_&1024)!=0&&CLEAN_FIXED) {
3964    for (int iRow=0;iRow<numberRows_;iRow++) {
3965      int iPivot=pivotVariable_[iRow];
3966      if (!flagged(iPivot)&&pivoted(iPivot)) {
3967        // carry on
3968        numberPrimalInfeasibilities_=-1;
3969        sumOfRelaxedPrimalInfeasibilities_ = 1.0;
3970        sumPrimalInfeasibilities_ = 1.0;
3971        break;
3972      }
3973    }
3974  }
3975  /* If we are primal feasible and any dual infeasibilities are on
3976     free variables then it is better to go to primal */
3977  if (!numberPrimalInfeasibilities_&&!numberDualInfeasibilitiesWithoutFree_&&
3978      numberDualInfeasibilities_)
3979    problemStatus_=10;
3980  // dual bound coming in
3981  //double saveDualBound = dualBound_;
3982  while (problemStatus_<=-3) {
3983    int cleanDuals=0;
3984    if (situationChanged!=0)
3985      cleanDuals=1;
3986    int numberChangedBounds=0;
3987    int doOriginalTolerance=0;
3988    if ( lastCleaned==numberIterations_)
3989      doOriginalTolerance=1;
3990    // check optimal
3991    // give code benefit of doubt
3992    if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
3993        sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
3994      // say optimal (with these bounds etc)
3995      numberDualInfeasibilities_ = 0;
3996      sumDualInfeasibilities_ = 0.0;
3997      numberPrimalInfeasibilities_ = 0;
3998      sumPrimalInfeasibilities_ = 0.0;
3999    }
4000    //if (dualFeasible()||problemStatus_==-4||(primalFeasible()&&!numberDualInfeasibilitiesWithoutFree_)) {
4001    if (dualFeasible()||problemStatus_==-4) {
4002      progress_.modifyObjective(objectiveValue_
4003                               -sumDualInfeasibilities_*dualBound_);
4004      // see if cutoff reached
4005      double limit = 0.0;
4006      getDblParam(ClpDualObjectiveLimit, limit);
4007      if(fabs(limit)<1.0e30&&objectiveValue()*optimizationDirection_>
4008         limit+1.0e-7+1.0e-8*fabs(limit)&&!numberAtFakeBound()) {
4009        //looks infeasible on objective
4010        if (perturbation_==101) {
4011          perturbation_=102; // stop any perturbations
4012          cleanDuals=1;
4013          // make sure fake bounds are back
4014          changeBounds(1,NULL,changeCost);
4015          createRim4(false);
4016          // make sure duals are current
4017          computeDuals(givenDuals);
4018          checkDualSolution();
4019        }
4020      }
4021      if (primalFeasible()&&!givenDuals) {
4022        // may be optimal - or may be bounds are wrong
4023        handler_->message(CLP_DUAL_BOUNDS,messages_)
4024          <<dualBound_
4025          <<CoinMessageEol;
4026        // save solution in case unbounded
4027        double * saveColumnSolution = NULL;
4028        double * saveRowSolution = NULL;
4029        bool inCbc = (specialOptions_&(0x01000000|16384))!=0;
4030        if (!inCbc) {
4031          saveColumnSolution = CoinCopyOfArray(columnActivityWork_,numberColumns_);
4032          saveRowSolution = CoinCopyOfArray(rowActivityWork_,numberRows_);
4033        }
4034        numberChangedBounds=changeBounds(0,rowArray_[3],changeCost);
4035        if (numberChangedBounds<=0&&!numberDualInfeasibilities_) {
4036          //looks optimal - do we need to reset tolerance
4037          if (perturbation_==101) {
4038            perturbation_=102; // stop any perturbations
4039            cleanDuals=1;
4040            // make sure fake bounds are back
4041            //computeObjectiveValue();
4042            changeBounds(1,NULL,changeCost);
4043            //computeObjectiveValue();
4044            createRim4(false);
4045            // make sure duals are current
4046            computeDuals(givenDuals);
4047            checkDualSolution();
4048#define DUAL_TRY_FASTER
4049#ifdef DUAL_TRY_FASTER
4050            if (numberDualInfeasibilities_) {
4051#endif
4052              numberChanged_=1; // force something to happen
4053#ifdef DUAL_TRY_FASTER
4054            } else {
4055              //double value = objectiveValue_;
4056              computeObjectiveValue(true);
4057              //printf("old %g new %g\n",value,objectiveValue_);
4058              //numberChanged_=1;
4059            }
4060#endif
4061          }
4062          if (lastCleaned<numberIterations_&&numberTimesOptimal_<4&&
4063              (numberChanged_||(specialOptions_&4096)==0)) {
4064            doOriginalTolerance=2;
4065            numberTimesOptimal_++;
4066            changeMade_++; // say something changed
4067            if (numberTimesOptimal_==1) {
4068              dualTolerance_ = dblParam_[ClpDualTolerance];
4069            } else {
4070              if (numberTimesOptimal_==2) {
4071                // better to have small tolerance even if slower
4072                factorization_->zeroTolerance(1.0e-15);
4073              }
4074              dualTolerance_ = dblParam_[ClpDualTolerance];
4075              dualTolerance_ *= pow(2.0,numberTimesOptimal_-1);
4076            }
4077            cleanDuals=2; // If nothing changed optimal else primal
4078          } else {
4079            problemStatus_=0; // optimal
4080            if (lastCleaned<numberIterations_&&numberChanged_) {
4081              handler_->message(CLP_SIMPLEX_GIVINGUP,messages_)
4082                <<CoinMessageEol;
4083            }
4084          }
4085        } else {
4086          cleanDuals=1;
4087          if (doOriginalTolerance==1) {
4088            // check unbounded
4089            // find a variable with bad dj
4090            int iSequence;
4091            int iChosen=-1;
4092            if (!inCbc) {
4093              double largest = 100.0*primalTolerance_;
4094              for (iSequence=0;iSequence<numberRows_+numberColumns_;
4095                   iSequence++) {
4096                double djValue = dj_[iSequence];
4097                double originalLo = originalLower(iSequence);
4098                double originalUp = originalUpper(iSequence);
4099                if (fabs(djValue)>fabs(largest)) {
4100                  if (getStatus(iSequence)!=basic) {
4101                    if (djValue>0&&originalLo<-1.0e20) {
4102                      if (djValue>fabs(largest)) {
4103                        largest=djValue;
4104                        iChosen=iSequence;
4105                      }
4106                    } else if (djValue<0&&originalUp>1.0e20) {
4107                      if (-djValue>fabs(largest)) {
4108                        largest=djValue;
4109                        iChosen=iSequence;
4110                      }
4111                    }
4112                  }
4113                }
4114              }
4115            }
4116            if (iChosen>=0) {
4117              int iSave=sequenceIn_;
4118              sequenceIn_=iChosen;
4119              unpack(rowArray_[1]);
4120              sequenceIn_ = iSave;
4121              // if dual infeasibilities then must be free vector so add in dual
4122              if (numberDualInfeasibilities_) {
4123                if (fabs(changeCost)>1.0e-5)
4124                  printf("Odd free/unbounded combo\n");
4125                changeCost += cost_[iChosen];
4126              }
4127              problemStatus_ = checkUnbounded(rowArray_[1],rowArray_[0],
4128                                              changeCost);
4129              rowArray_[1]->clear();
4130            } else {
4131              problemStatus_=-3;
4132            }
4133            if (problemStatus_==2&&perturbation_==101) {
4134              perturbation_=102; // stop any perturbations
4135              cleanDuals=1;
4136              createRim4(false);
4137              problemStatus_=-1;
4138            }
4139            if (problemStatus_==2) {
4140              // it is unbounded - restore solution
4141              // but first add in changes to non-basic
4142              int iColumn;
4143              double * original = columnArray_[0]->denseVector();
4144              for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4145                if(getColumnStatus(iColumn)!= basic)
4146                  ray_[iColumn] +=
4147                    saveColumnSolution[iColumn]-original[iColumn];
4148                columnActivityWork_[iColumn] = original[iColumn];
4149              }
4150              CoinMemcpyN(saveRowSolution,numberRows_,
4151                                rowActivityWork_);
4152            }
4153          } else {
4154            doOriginalTolerance=2;
4155            rowArray_[0]->clear();
4156          }
4157        }
4158        delete [] saveColumnSolution;
4159        delete [] saveRowSolution;
4160      }
4161      if (problemStatus_==-4||problemStatus_==-5) {
4162        // may be infeasible - or may be bounds are wrong
4163        numberChangedBounds=changeBounds(0,NULL,changeCost);
4164        /* Should this be here as makes no difference to being feasible.
4165           But seems to make a difference to run times. */
4166        if (perturbation_==101&&0) {
4167          perturbation_=102; // stop any perturbations
4168          cleanDuals=1;
4169          numberChangedBounds=1;
4170          // make sure fake bounds are back
4171          changeBounds(1,NULL,changeCost);
4172          createRim4(false);
4173        }
4174        if (numberChangedBounds<=0||dualBound_>1.0e20||
4175            (largestPrimalError_>1.0&&dualBound_>1.0e17)) {
4176          problemStatus_=1; // infeasible
4177          if (perturbation_==101) {
4178            perturbation_=102; // stop any perturbations
4179            //cleanDuals=1;
4180            //numberChangedBounds=1;
4181            //createRim4(false);
4182          }
4183        } else {
4184          problemStatus_=-1; //iterate
4185          cleanDuals=1;
4186          if (numberChangedBounds<=0)
4187            doOriginalTolerance=2;
4188          // and delete ray which has been created
4189          delete [] ray_;
4190          ray_ = NULL;
4191        }
4192
4193      }
4194    } else {
4195      cleanDuals=1;
4196    }
4197    if (problemStatus_<0) {
4198      if (doOriginalTolerance==2) {
4199        // put back original tolerance
4200        lastCleaned=numberIterations_;
4201        numberChanged_ =0; // Number of variables with changed costs
4202        handler_->message(CLP_DUAL_ORIGINAL,messages_)
4203          <<CoinMessageEol;
4204        perturbation_=102; // stop any perturbations
4205#if 0
4206        double * xcost = new double[numberRows_+numberColumns_];
4207        double * xlower = new double[numberRows_+numberColumns_];
4208        double * xupper = new double[numberRows_+numberColumns_];
4209        double * xdj = new double[numberRows_+numberColumns_];
4210        double * xsolution = new double[numberRows_+numberColumns_];
4211 CoinMemcpyN(cost_,(numberRows_+numberColumns_),xcost);
4212 CoinMemcpyN(lower_,(numberRows_+numberColumns_),xlower);
4213 CoinMemcpyN(upper_,(numberRows_+numberColumns_),xupper);
4214 CoinMemcpyN(dj_,(numberRows_+numberColumns_),xdj);
4215 CoinMemcpyN(solution_,(numberRows_+numberColumns_),xsolution);
4216#endif
4217        createRim4(false);
4218        // make sure duals are current
4219        computeDuals(givenDuals);
4220        checkDualSolution();
4221#if 0
4222        int i;
4223        for (i=0;i<numberRows_+numberColumns_;i++) {
4224          if (cost_[i]!=xcost[i])
4225            printf("** %d old cost %g new %g sol %g\n",
4226                   i,xcost[i],cost_[i],solution_[i]);
4227          if (lower_[i]!=xlower[i])
4228            printf("** %d old lower %g new %g sol %g\n",
4229                   i,xlower[i],lower_[i],solution_[i]);
4230          if (upper_[i]!=xupper[i])
4231            printf("** %d old upper %g new %g sol %g\n",
4232                   i,xupper[i],upper_[i],solution_[i]);
4233          if (dj_[i]!=xdj[i])
4234            printf("** %d old dj %g new %g sol %g\n",
4235                   i,xdj[i],dj_[i],solution_[i]);
4236          if (solution_[i]!=xsolution[i])
4237            printf("** %d old solution %g new %g sol %g\n",
4238                   i,xsolution[i],solution_[i],solution_[i]);
4239        }
4240        //delete [] xcost;
4241        //delete [] xupper;
4242        //delete [] xlower;
4243        //delete [] xdj;
4244        //delete [] xsolution;
4245#endif
4246        // put back bounds as they were if was optimal
4247        if (doOriginalTolerance==2&&cleanDuals!=2) {
4248          changeMade_++; // say something changed
4249          /* We may have already changed some bounds in this function
4250             so save numberFake_ and add in.
4251
4252             Worst that can happen is that we waste a bit of time  - but it must be finite.
4253          */
4254          int saveNumberFake = numberFake_;
4255          changeBounds(1,NULL,changeCost);
4256          numberFake_ += saveNumberFake;
4257          cleanDuals=2;
4258          //cleanDuals=1;
4259        }
4260#if 0
4261        //int i;
4262        for (i=0;i<numberRows_+numberColumns_;i++) {
4263          if (cost_[i]!=xcost[i])
4264            printf("** %d old cost %g new %g sol %g\n",
4265                   i,xcost[i],cost_[i],solution_[i]);
4266          if (lower_[i]!=xlower[i])
4267            printf("** %d old lower %g new %g sol %g\n",
4268                   i,xlower[i],lower_[i],solution_[i]);
4269          if (upper_[i]!=xupper[i])
4270            printf("** %d old upper %g new %g sol %g\n",
4271                   i,xupper[i],upper_[i],solution_[i]);
4272          if (dj_[i]!=xdj[i])
4273            printf("** %d old dj %g new %g sol %g\n",
4274                   i,xdj[i],dj_[i],solution_[i]);
4275          if (solution_[i]!=xsolution[i])
4276            printf("** %d old solution %g new %g sol %g\n",
4277                   i,xsolution[i],solution_[i],solution_[i]);
4278        }
4279        delete [] xcost;
4280        delete [] xupper;
4281        delete [] xlower;
4282        delete [] xdj;
4283        delete [] xsolution;
4284#endif
4285      }
4286      if (cleanDuals==1||(cleanDuals==2&&!numberDualInfeasibilities_)) {
4287        // make sure dual feasible
4288        // look at all rows and columns
4289        rowArray_[0]->clear();
4290        columnArray_[0]->clear();
4291        double objectiveChange=0.0;
4292#if 0
4293        double * xcost = new double[numberRows_+numberColumns_];
4294        double * xlower = new double[numberRows_+numberColumns_];
4295        double * xupper = new double[numberRows_+numberColumns_];
4296        double * xdj = new double[numberRows_+numberColumns_];
4297        double * xsolution = new double[numberRows_+numberColumns_];
4298 CoinMemcpyN(cost_,(numberRows_+numberColumns_),xcost);
4299 CoinMemcpyN(lower_,(numberRows_+numberColumns_),xlower);
4300 CoinMemcpyN(upper_,(numberRows_+numberColumns_),xupper);
4301 CoinMemcpyN(dj_,(numberRows_+numberColumns_),xdj);
4302 CoinMemcpyN(solution_,(numberRows_+numberColumns_),xsolution);
4303#endif
4304        if (givenDuals)
4305          dualTolerance_=1.0e50;
4306        updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[1],
4307          0.0,objectiveChange,true);
4308        dualTolerance_=saveTolerance;
4309#if 0
4310        int i;
4311        for (i=0;i<numberRows_+numberColumns_;i++) {
4312          if (cost_[i]!=xcost[i])
4313            printf("** %d old cost %g new %g sol %g\n",
4314                   i,xcost[i],cost_[i],solution_[i]);
4315          if (lower_[i]!=xlower[i])
4316            printf("** %d old lower %g new %g sol %g\n",
4317                   i,xlower[i],lower_[i],solution_[i]);
4318          if (upper_[i]!=xupper[i])
4319            printf("** %d old upper %g new %g sol %g\n",
4320                   i,xupper[i],upper_[i],solution_[i]);
4321          if (dj_[i]!=xdj[i])
4322            printf("** %d old dj %g new %g sol %g\n",
4323                   i,xdj[i],dj_[i],solution_[i]);
4324          if (solution_[i]!=xsolution[i])
4325            printf("** %d old solution %g new %g sol %g\n",
4326                   i,xsolution[i],solution_[i],solution_[i]);
4327        }
4328        delete [] xcost;
4329        delete [] xupper;
4330        delete [] xlower;
4331        delete [] xdj;
4332        delete [] xsolution;
4333#endif
4334        // for now - recompute all
4335        gutsOfSolution(NULL,NULL);
4336        if (givenDuals)
4337          dualTolerance_=1.0e50;
4338        updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[1],
4339          0.0,objectiveChange,true);
4340        dualTolerance_=saveTolerance;
4341        //assert(numberDualInfeasibilitiesWithoutFree_==0);
4342        if (numberDualInfeasibilities_) {
4343          if (numberPrimalInfeasibilities_||numberPivots) {
4344            problemStatus_=-1; // carry on as normal
4345          } else {
4346            problemStatus_=10; // try primal
4347#if COIN_DEVELOP>1
4348            printf("returning at %d\n",__LINE__);
4349#endif
4350          }
4351        } else if (situationChanged==2) {
4352          problemStatus_=-1; // carry on as normal
4353        }
4354        situationChanged=0;
4355      } else {
4356        // iterate
4357        if (cleanDuals!=2) {
4358          problemStatus_=-1;
4359        } else {
4360          problemStatus_ = 10; // try primal
4361#if COIN_DEVELOP>2
4362          printf("returning at %d\n",__LINE__);
4363#endif
4364        }
4365      }
4366    }
4367  }
4368  // unflag all variables (we may want to wait a bit?)
4369  if ((tentativeStatus!=-2&&tentativeStatus!=-1)&&unflagVariables) {
4370    int iRow;
4371    int numberFlagged=0;
4372    for (iRow=0;iRow<numberRows_;iRow++) {
4373      int iPivot=pivotVariable_[iRow];
4374      if (flagged(iPivot)) {
4375        numberFlagged++;
4376        clearFlagged(iPivot);
4377      }
4378    }
4379#ifdef COIN_DEVELOP
4380    if (numberFlagged) {
4381      printf("unflagging %d variables - tentativeStatus %d probStat %d ninf %d nopt %d\n",numberFlagged,tentativeStatus,
4382             problemStatus_,numberPrimalInfeasibilities_,
4383             numberTimesOptimal_);
4384    }
4385#endif
4386    unflagVariables = numberFlagged>0;
4387    if (numberFlagged&&!numberPivots) {
4388      /* looks like trouble as we have not done any iterations.
4389         Try changing pivot tolerance then give it a few goes and give up */
4390      if (factorization_->pivotTolerance()<0.9) {
4391        factorization_->pivotTolerance(0.99);
4392        problemStatus_=-1;
4393      } else if (numberTimesOptimal_<3) {
4394        numberTimesOptimal_++;
4395        problemStatus_=-1;
4396      } else {
4397        unflagVariables=false;
4398        //secondaryStatus_ = 1; // and say probably infeasible
4399        // try primal
4400        problemStatus_=10;
4401#if COIN_DEVELOP>1
4402        printf("returning at %d\n",__LINE__);
4403#endif
4404      }
4405    }
4406  }
4407  if (problemStatus_<0) {
4408    if (type==0||type==1) {
4409      if (!type) {
4410        // create save arrays
4411        delete [] saveStatus_;
4412        delete [] savedSolution_;
4413        saveStatus_ = new unsigned char [numberRows_+numberColumns_];
4414        savedSolution_ = new double [numberRows_+numberColumns_];
4415      }
4416      // save arrays
4417      CoinMemcpyN(status_,numberColumns_+numberRows_,saveStatus_);
4418      CoinMemcpyN(rowActivityWork_,
4419                  numberRows_,savedSolution_+numberColumns_);
4420      CoinMemcpyN(columnActivityWork_,numberColumns_,savedSolution_);
4421      // save extra stuff
4422      int dummy;
4423      matrix_->generalExpanded(this,5,dummy);
4424    }
4425    if (weightsSaved) {
4426      // restore weights (if saved) - also recompute infeasibility list
4427      if (tentativeStatus>-3)
4428        dualRowPivot_->saveWeights(this,(type <2) ? 2 : 4);
4429      else
4430        dualRowPivot_->saveWeights(this,3);
4431    }
4432  }
4433  // see if cutoff reached
4434  double limit = 0.0;
4435  getDblParam(ClpDualObjectiveLimit, limit);
4436#if 0
4437  if(fabs(limit)<1.0e30&&objectiveValue()*optimizationDirection_>
4438     limit+100.0) {
4439    printf("lim %g obj %g %g - wo perturb %g sum dual %g\n",
4440           limit,objectiveValue_,objectiveValue(),computeInternalObjectiveValue(),sumDualInfeasibilities_);
4441  }
4442#endif
4443  if(fabs(limit)<1.0e30&&objectiveValue()*optimizationDirection_>
4444     limit&&!numberAtFakeBound()) {
4445    bool looksInfeasible = !numberDualInfeasibilities_;
4446    if (objectiveValue()*optimizationDirection_>limit+fabs(0.1*limit)+1.0e2*sumDualInfeasibilities_+1.0e4&&
4447        sumDualInfeasibilities_<largestDualError_&&numberIterations_>0.5*numberRows_+1000)
4448      looksInfeasible=true;
4449    if (looksInfeasible) {
4450      // Even if not perturbed internal costs may have changed
4451      // be careful
4452      if (numberIterations_) {
4453        if(computeInternalObjectiveValue()>limit) {
4454          problemStatus_=1;
4455          secondaryStatus_ = 1; // and say was on cutoff
4456        }
4457      } else {
4458        problemStatus_=1;
4459        secondaryStatus_ = 1; // and say was on cutoff
4460      }
4461    }
4462  }
4463  // If we are in trouble and in branch and bound give up
4464  if ((specialOptions_&1024)!=0) {
4466    if (largestPrimalError_*largestDualError_>1.0e2) {
4468    } else if (largestPrimalError_>1.0e-2
4469        &&objectiveValue_>CoinMin(1.0e15,1.0e3*limit)) {
4471    }
4473      if (factorization_->pivotTolerance()<0.9) {
4474        // up tolerance
4475        factorization_->pivotTolerance(CoinMin(factorization_->pivotTolerance()*1.05+0.02,0.91));
4476      } else if (numberIterations_>10000) {
4477        if (handler_->logLevel()>2)
4479        problemStatus_=1;
4480        secondaryStatus_ = 1; // and say was on cutoff
4481      } else if (largestPrimalError_>1.0e5) {
4482        {
4483          int iBigB=-1;
4484          double bigB=0.0;
4485          int iBigN=-1;
4486          double bigN=0.0;
4487          for (int i=0;i<numberRows_+numberColumns_;i++) {
4488            double value = fabs(solution_[i]);
4489            if (getStatus(i)==basic) {
4490              if (value>bigB) {
4491                bigB= value;
4492                iBigB=i;
4493              }
4494            } else {
4495              if (value>bigN) {
4496                bigN= value;
4497                iBigN=i;
4498              }
4499            }
4500          }
4501#ifdef CLP_INVESTIGATE
4502          if (bigB>1.0e8||bigN>1.0e8) {
4503            if (handler_->logLevel()>0)
4504              printf("it %d - basic %d %g, nonbasic %d %g\n",
4505                     numberIterations_,iBigB,bigB,iBigN,bigN);
4506          }
4507#endif
4508        }
4509#if COIN_DEVELOP!=2
4510        if (handler_->logLevel()>2)
4511#endif
4513        allSlackBasis(true);
4514        problemStatus_=10;
4515      }
4516    }
4517  }
4519    problemStatus_=4; // unknown
4520  }
4521  lastGoodIteration_ = numberIterations_;
4523    moreSpecialOptions_ &= ~16; // clear check accuracy flag
4524  if (problemStatus_<0) {
4525    sumDualInfeasibilities_=realDualInfeasibilities; // back to say be careful
4526    if (sumDualInfeasibilities_)
4527      numberDualInfeasibilities_=1;
4528  }
4529#if 1
4530  double thisObj = progress_.lastObjective(0);
4531  double lastObj = progress_.lastObjective(1);
4532  if (lastObj>thisObj+1.0e-4*CoinMax(fabs(thisObj),fabs(lastObj))+1.0e-4
4533      &&givenDuals==NULL&&firstFree_<0) {
4534    int maxFactor = factorization_->maximumPivots();
4535    if (maxFactor>10) {
4536      if (forceFactorization_<0)
4537        forceFactorization_= maxFactor;
4538      forceFactorization_ = CoinMax(1,(forceFactorization_>>1));
4539      //printf("Reducing factorization frequency\n");
4540    }
4541  }
4542#endif
4543  // Allow matrices to be sorted etc
4544  int fake=-999; // signal sort
4545  matrix_->correctSequence(this,fake,fake);
4546  if (alphaAccuracy_>0.0)
4547      alphaAccuracy_=1.0;
4548}
4549/* While updateDualsInDual sees what effect is of flip
4550   this does actual flipping.
4551   If change >0.0 then value in array >0.0 => from lower to upper
4552*/
4553void
4554ClpSimplexDual::flipBounds(CoinIndexedVector * rowArray,
4555                  CoinIndexedVector * columnArray,
4556                  double change)
4557{
4558  int number;
4559  int * which;
4560
4561  int iSection;
4562
4563  for (iSection=0;iSection<2;iSection++) {
4564    int i;
4565    double * solution = solutionRegion(iSection);
4566    double * lower = lowerRegion(iSection);
4567    double * upper = upperRegion(iSection);
4569    if (!iSection) {
4570      number = rowArray->getNumElements();
4571      which = rowArray->getIndices();
4573    } else {
4574      number = columnArray->getNumElements();
4575      which = columnArray->getIndices();
4577    }
4578
4579    for (i=0;i<number;i++) {
4580      int iSequence = which[i];
4582
4583      switch(status) {
4584
4585      case basic:
4586      case isFree:
4587      case superBasic:
4588      case ClpSimplex::isFixed:
4589        break;
4590      case atUpperBound:
4591        // to lower bound
4593        solution[iSequence] = lower[iSequence];
4594        break;
4595      case atLowerBound:
4596        // to upper bound
4598        solution[iSequence] = upper[iSequence];
4599        break;
4600      }
4601    }
4602  }
4603  rowArray->setNumElements(0);
4604  columnArray->setNumElements(0);
4605}
4606// Restores bound to original bound
4607void
4608ClpSimplexDual::originalBound( int iSequence)
4609{
4610  if (getFakeBound(iSequence)!=noFake) {
4611    numberFake_--;;
4612    setFakeBound(iSequence,noFake);
4613    if (auxiliaryModel_) {
4614      // just copy back
4615      lower_[iSequence]=auxiliaryModel_->lowerRegion()[iSequence+numberRows_+numberColumns_];
4616      upper_[iSequence]=auxiliaryModel_->upperRegion()[iSequence+numberRows_+numberColumns_];
4617      return;
4618    }
4619    if (iSequence>=numberColumns_) {
4620      // rows
4621      int iRow = iSequence-numberColumns_;
4622      rowLowerWork_[iRow]=rowLower_[iRow];
4623      rowUpperWork_[iRow]=rowUpper_[iRow];
4624      if (rowScale_) {
4625        if (rowLowerWork_[iRow]>-1.0e50)
4626          rowLowerWork_[iRow] *= rowScale_[iRow]*rhsScale_;
4627        if (rowUpperWork_[iRow]<1.0e50)
4628          rowUpperWork_[iRow] *= rowScale_[iRow]*rhsScale_;
4629      } else if (rhsScale_!=1.0) {
4630        if (rowLowerWork_[iRow]>-1.0e50)
4631          rowLowerWork_[iRow] *= rhsScale_;
4632        if (rowUpperWork_[iRow]<1.0e50)
4633          rowUpperWork_[iRow] *= rhsScale_;
4634      }
4635    } else {
4636      // columns
4637      columnLowerWork_[iSequence]=columnLower_[iSequence];
4638      columnUpperWork_[iSequence]=columnUpper_[iSequence];
4639      if (rowScale_) {
4640        double multiplier = 1.0/columnScale_[iSequence];
4641        if (columnLowerWork_[iSequence]>-1.0e50)
4642          columnLowerWork_[iSequence] *= multiplier*rhsScale_;
4643        if (columnUpperWork_[iSequence]<1.0e50)
4644          columnUpperWork_[iSequence] *= multiplier*rhsScale_;
4645      } else if (rhsScale_!=1.0) {
4646        if (columnLowerWork_[iSequence]>-1.0e50)
4647          columnLowerWork_[iSequence] *= rhsScale_;
4648        if (columnUpperWork_[iSequence]<1.0e50)
4649          columnUpperWork_[iSequence] *= rhsScale_;
4650      }
4651    }
4652  }
4653}
4654/* As changeBounds but just changes new bounds for a single variable.
4655   Returns true if change */
4656bool
4657ClpSimplexDual::changeBound( int iSequence)
4658{
4659  // old values
4660  double oldLower=lower_[iSequence];
4661  double oldUpper=upper_[iSequence];
4662  double value=solution_[iSequence];
4663  bool modified=false;
4664  originalBound(iSequence);
4665  // original values
4666  double lowerValue=lower_[iSequence];
4667  double upperValue=upper_[iSequence];
4668  // back to altered values
4669  lower_[iSequence] = oldLower;
4670  upper_[iSequence] = oldUpper;
4671  assert (getFakeBound(iSequence)==noFake);
4672  //if (getFakeBound(iSequence)!=noFake)
4673  //numberFake_--;;
4674  if (value==oldLower) {
4675    if (upperValue > oldLower + dualBound_) {
4676      upper_[iSequence]=oldLower+dualBound_;
4677      setFakeBound(iSequence,upperFake);
4678      modified=true;
4679      numberFake_++;
4680    }
4681  } else if (value==oldUpper) {
4682    if (lowerValue < oldUpper - dualBound_) {
4683      lower_[iSequence]=oldUpper-dualBound_;
4684      setFakeBound(iSequence,lowerFake);
4685      modified=true;
4686      numberFake_++;
4687    }
4688  } else {
4689    assert(value==oldLower||value==oldUpper);
4690  }
4691  return modified;
4692}
4693// Perturbs problem
4694int
4695ClpSimplexDual::perturb()
4696{
4697  if (perturbation_>100)
4699  if (perturbation_==100)
4700    perturbation_=50; // treat as normal
4701  int savePerturbation = perturbation_;
4702  bool modifyRowCosts=false;
4703  // dual perturbation
4704  double perturbation=1.0e-20;
4705  // maximum fraction of cost to perturb
4706  double maximumFraction = 1.0e-5;
4707  double constantPerturbation = 100.0*dualTolerance_;
4708  int maxLength=0;
4709  int minLength=numberRows_;
4710  double averageCost = 0.0;
4711#if 0
4712  // look at element range
4713  double smallestNegative;
4714  double largestNegative;
4715  double smallestPositive;
4716  double largestPositive;
4717  matrix_->rangeOfElements(smallestNegative, largestNegative,
4718                           smallestPositive, largestPositive);
4719  smallestPositive = CoinMin(fabs(smallestNegative),smallestPositive);
4720  largestPositive = CoinMax(fabs(largestNegative),largestPositive);
4721  double elementRatio = largestPositive/smallestPositive;
4722#endif
4723  int numberNonZero=0;
4724  if (!numberIterations_&&perturbation_>=50) {
4725    // See if we need to perturb
4726    double * sort = new double[numberColumns_];
4727    // Use objective BEFORE scaling
4728    const double * obj = objective();
4729    int i;
4730    for (i=0;i<numberColumns_;i++) {
4731      double value = fabs(obj[i]);
4732      sort[i]=value;
4733      averageCost += value;
4734      if (value)
4735        numberNonZero++;
4736    }
4737    if (numberNonZero)
4738      averageCost /= (double) numberNonZero;
4739    else
4740      averageCost = 1.0;
4741    std::sort(sort,sort+numberColumns_);
4742    int number=1;
4743    double last = sort[0];
4744    for (i=1;i<numberColumns_;i++) {
4745      if (last!=sort[i])
4746        number++;
4747      last=sort[i];
4748    }
4749    delete [] sort;
4750    if (!numberNonZero)
4751      return 1; // safer to use primal
4752#if 0
4753    printf("nnz %d percent %d",number,(number*100)/numberColumns_);
4754    if (number*4>numberColumns_)
4755      printf(" - Would not perturb\n");
4756    else
4757      printf(" - Would perturb\n");
4758    //exit(0);
4759#endif
4760    //printf("ratio number diff costs %g, element ratio %g\n",((double)number)/((double) numberColumns_),
4761    //                                                                elementRatio);
4762    //number=0;
4763    //if (number*4>numberColumns_||elementRatio>1.0e12) {
4764    if (number*4>numberColumns_) {
4765      perturbation_=100;
4766      return 0; // good enough
4767    }
4768  }
4769  int iColumn;
4770  const int * columnLength = matrix_->getVectorLengths();
4771  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4772    if (columnLowerWork_[iColumn]<columnUpperWork_[iColumn]) {
4773      int length = columnLength[iColumn];
4774      if (length>2) {
4775        maxLength = CoinMax(maxLength,length);
4776        minLength = CoinMin(minLength,length);
4777      }
4778    }
4779  }
4780  // If > 70 then do rows
4781  if (perturbation_>=70) {
4782    modifyRowCosts=true;
4783    perturbation_ -= 20;
4784    printf("Row costs modified, ");
4785  }
4786  bool uniformChange=false;
4787  bool inCbcOrOther = (specialOptions_&0x03000000)!=0;
4788  if (perturbation_>50) {
4789    // Experiment
4790    // maximumFraction could be 1.0e-10 to 1.0
4791    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};
4792    int whichOne = perturbation_-51;
4793    //if (inCbcOrOther&&whichOne>0)
4794    //whichOne--;
4795    maximumFraction = m[CoinMin(whichOne,10)];
4796  } else if (inCbcOrOther) {
4797    //maximumFraction = 1.0e-6;
4798  }
4799  int iRow;
4800  double smallestNonZero=1.0e100;
4801  numberNonZero=0;
4802  if (perturbation_>=50) {
4803    perturbation = 1.0e-8;
4804    bool allSame=true;
4805    double lastValue=0.0;
4806    for (iRow=0;iRow<numberRows_;iRow++) {
4807      double lo = rowLowerWork_[iRow];
4808      double up = rowUpperWork_[iRow];
4809      if (lo<up) {
4810        double value = fabs(rowObjectiveWork_[iRow]);
4811        perturbation = CoinMax(perturbation,value);
4812        if (value) {
4813          modifyRowCosts=true;
4814          smallestNonZero = CoinMin(smallestNonZero,value);
4815        }
4816      }
4817      if (lo&&lo>-1.0e10) {
4818        numberNonZero++;
4819        lo=fabs(lo);
4820        if (!lastValue)
4821          lastValue=lo;
4822        else if (fabs(lo-lastValue)>1.0e-7)
4823          allSame=false;
4824      }
4825      if (up&&up<1.0e10) {
4826        numberNonZero++;
4827        up=fabs(up);
4828        if (!lastValue)
4829          lastValue=up;
4830        else if (fabs(up-lastValue)>1.0e-7)
4831          allSame=false;
4832      }
4833    }
4834    double lastValue2=0.0;
4835    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4836      double lo = columnLowerWork_[iColumn];
4837      double up = columnUpperWork_[iColumn];
4838      if (lo<up) {
4839        double value =
4840          fabs(objectiveWork_[iColumn]);
4841        perturbation = CoinMax(perturbation,value);
4842        if (value) {
4843          smallestNonZero = CoinMin(smallestNonZero,value);
4844        }
4845      }
4846      if (lo&&lo>-1.0e10) {
4847        //numberNonZero++;
4848        lo=fabs(lo);
4849        if (!lastValue2)
4850          lastValue2=lo;
4851        else if (fabs(lo-lastValue2)>1.0e-7)
4852          allSame=false;
4853      }
4854      if (up&&up<1.0e10) {
4855        //numberNonZero++;
4856        up=fabs(up);
4857        if (!lastValue2)
4858          lastValue2=up;
4859        else if (fabs(up-lastValue2)>1.0e-7)
4860          allSame=false;
4861      }
4862    }
4863    if (allSame) {
4864      // Check elements
4865      double smallestNegative;
4866      double largestNegative;
4867      double smallestPositive;
4868      double largestPositive;
4869      matrix_->rangeOfElements(smallestNegative,largestNegative,
4870                               smallestPositive,largestPositive);
4871      if (smallestNegative==largestNegative&&
4872          smallestPositive==largestPositive) {
4873        // Really hit perturbation
4876      }
4877    }
4878    perturbation = CoinMin(perturbation,smallestNonZero/maximumFraction);
4879  } else {
4880    // user is in charge
4881    maximumFraction = 1.0e-1;
4882    // but some experiments
4883    if (perturbation_<=-900) {
4884      modifyRowCosts=true;
4885      perturbation_ += 1000;
4886      printf("Row costs modified, ");
4887    }
4888    if (perturbation_<=-10) {
4889      perturbation_ += 10;
4890      maximumFraction = 1.0;
4891      if ((-perturbation_)%100>=10) {
4892        uniformChange=true;
4893        perturbation_+=20;
4894      }
4895      while (perturbation_<-10) {
4896        perturbation_ += 100;
4897        maximumFraction *= 1.0e-1;
4898      }
4899    }
4900    perturbation = pow(10.0,perturbation_);
4901  }
4902  double largestZero=0.0;
4903  double largest=0.0;
4904  double largestPerCent=0.0;
4905  // modify costs
4906  bool printOut=(handler_->logLevel()==63);
4907  printOut=false;
4908  if (modifyRowCosts) {
4909    for (iRow=0;iRow<numberRows_;iRow++) {
4910      if (rowLowerWork_[iRow]<rowUpperWork_[iRow]) {
4911        double value = perturbation;
4912        double currentValue = rowObjectiveWork_[iRow];
4913        value = CoinMin(value,maximumFraction*(fabs(currentValue)+1.0e-1*perturbation+1.0e-3));
4914        if (rowLowerWork_[iRow]>-largeValue_) {
4915          if (fabs(rowLowerWork_[iRow])<fabs(rowUpperWork_[iRow]))
4916            value *= randomNumberGenerator_.randomDouble();
4917          else
4918            value *= -randomNumberGenerator_.randomDouble();
4919        } else if (rowUpperWork_[iRow]<largeValue_) {
4920          value *= -randomNumberGenerator_.randomDouble();
4921        } else {
4922          value=0.0;
4923        }
4924        if (currentValue) {
4925          largest = CoinMax(largest,fabs(value));
4926          if (fabs(value)>fabs(currentValue)*largestPerCent)
4927            largestPerCent=fabs(value/currentValue);
4928        } else {
4929          largestZero=CoinMax(largestZero,fabs(value));
4930        }
4931        if (printOut)
4932          printf("row %d cost %g change %g\n",iRow,rowObjectiveWork_[iRow],value);
4933        rowObjectiveWork_[iRow] += value;
4934      }
4935    }
4936  }
4937  // 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};
4938  // 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};
4939  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};
4940  //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};
4941  double extraWeight=10.0;
4942  // Scale back if wanted
4943  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};
4944  if (constantPerturbation<99.0*dualTolerance_) {
4945    perturbation *= 0.1;
4946    extraWeight=0.5;
4947    memcpy(weight,weight2,sizeof(weight2));
4948  }
4949  // adjust weights if all columns long
4950  double factor=1.0;
4951  if (maxLength) {
4952    factor = 3.0/(double) minLength;
4953  }
4954  // Make variables with more elements more expensive
4955  const double m1 = 0.5;
4956  double smallestAllowed = CoinMin(1.0e-2*dualTolerance_,maximumFraction);
4957  //double largestAllowed = CoinMax(1.0e3*dualTolerance_,maximumFraction*10.0*averageCost);
4958  double largestAllowed = CoinMax(1.0e3*dualTolerance_,maximumFraction*averageCost);
4959  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4960    if (columnLowerWork_[iColumn]<columnUpperWork_[iColumn]&&getStatus(iColumn)!=basic) {
4961      double value = perturbation;
4962      double currentValue = objectiveWork_[iColumn];
4963      value = CoinMin(value,constantPerturbation+maximumFraction*(fabs(currentValue)+1.0e-1*perturbation+1.0e-8));
4964      //value = CoinMin(value,constantPerturbation;+maximumFraction*fabs(currentValue));
4965      double value2 = constantPerturbation+1.0e-1*smallestNonZero;
4966      if (uniformChange) {
4967        value = maximumFraction;
4968        value2=maximumFraction;
4969      }
4970      if (columnLowerWork_[iColumn]>-largeValue_) {
4971        if (fabs(columnLowerWork_[iColumn])<
4972            fabs(columnUpperWork_[iColumn])) {
4973          value *= (1.0-m1 +m1*randomNumberGenerator_.randomDouble());
4974          value2 *= (1.0-m1 +m1*randomNumberGenerator_.randomDouble());
4975        } else {
4976          //value *= -(1.0-m1+m1*randomNumberGenerator_.randomDouble());
4977          //value2 *= -(1.0-m1+m1*randomNumberGenerator_.randomDouble());
4978          value=0.0;
4979        }
4980      } else if (columnUpperWork_[iColumn]<largeValue_) {
4981        value *= -(1.0-m1+m1*randomNumberGenerator_.randomDouble());
4982        value2 *= -(1.0-m1+m1*randomNumberGenerator_.randomDouble());
4983      } else {
4984        value=0.0;
4985      }
4986      if (value) {
4987        int length = columnLength[iColumn];
4988        if (length>3) {
4989          length = (int) ((double) length * factor);
4990          length = CoinMax(3,length);
4991        }
4992        double multiplier;
4993#if 1
4994        if (length<10)
4995          multiplier=weight[length];
4996        else
4997          multiplier = weight[10];
4998#else
4999        if (length<10)
5000          multiplier=weight[length];
5001        else
5002          multiplier = weight[10]+extraWeight*(length-10);
5003        multiplier *= 0.5;
5004#endif
5005        value *= multiplier;
5006        value = CoinMin(value,value2);
5007        if (savePerturbation<50||savePerturbation>60) {
5008          if (fabs(value)<=dualTolerance_)
5009            value=0.0;
5010        } else if (value) {
5011          // get in range
5012          if (fabs(value)<=smallestAllowed) {
5013            value *= 10.0;
5014            while (fabs(value)<=smallestAllowed)
5015              value *= 10.0;
5016          } else if (fabs(value)>largestAllowed) {
5017            value *= 0.1;
5018            while (fabs(value)>largestAllowed)
5019              value *= 0.1;
5020          }
5021        }
5022        if (currentValue) {
5023          largest = CoinMax(largest,fabs(value));
5024          if (fabs(value)>fabs(currentValue)*largestPerCent)
5025            largestPerCent=fabs(value/currentValue);
5026        } else {
5027          largestZero=CoinMax(largestZero,fabs(value));
5028        }
5029        // but negative if at ub
5030        if (getStatus(iColumn)==atUpperBound)
5031          value = -value;
5032        if (printOut)
5033          printf("col %d cost %g change %g\n",iColumn,objectiveWork_[iColumn],value);
5034        objectiveWork_[iColumn] += value;
5035      }
5036    }
5037  }
5038  handler_->message(CLP_SIMPLEX_PERTURB,messages_)
5039    <<100.0*maximumFraction<<perturbation<<largest<<100.0*largestPerCent<<largestZero
5040    <<CoinMessageEol;
5041  // and zero changes
5042  //int nTotal = numberRows_+numberColumns_;
5043  //CoinZeroN(cost_+nTotal,nTotal);
5044  // say perturbed
5045  perturbation_=101;
5046  return 0;
5047}
5048/* For strong branching.  On input lower and upper are new bounds
5049   while on output they are change in objective function values
5050   (>1.0e50 infeasible).
5051   Return code is 0 if nothing interesting, -1 if infeasible both
5052   ways and +1 if infeasible one way (check values to see which one(s))
5053   Returns -2 if bad factorization
5054*/
5055int ClpSimplexDual::strongBranching(int numberVariables,const int * variables,
5056                                    double * newLower, double * newUpper,
5057                                    double ** outputSolution,
5058                                    int * outputStatus, int * outputIterations,
5059                                    bool stopOnFirstInfeasible,
5060                                    bool alwaysFinish,
5061                                    int startFinishOptions)
5062{
5063  int i;
5064  int returnCode=0;
5065  double saveObjectiveValue = objectiveValue_;
5066  algorithm_ = -1;
5067
5068  //scaling(false);
5069
5070  // put in standard form (and make row copy)
5071  // create modifiable copies of model rim and do optional scaling
5072  createRim(7+8+16+32,true,startFinishOptions);
5073
5074  // change newLower and newUpper if scaled
5075
5076  // Do initial factorization
5077  // and set certain stuff
5078  // We can either set increasing rows so ...IsBasic gives pivot row
5079  // or we can just increment iBasic one by one
5080  // for now let ...iBasic give pivot row
5081  int useFactorization=false;
5082  if ((startFinishOptions&2)!=0&&(whatsChanged_&(2+512))==2+512)
5083    useFactorization=true; // Keep factorization if possible
5084  // switch off factorization if bad
5085  if (pivotVariable_[0]<0)
5086    useFactorization=false;
5087  if (!useFactorization||factorization_->numberRows()!=numberRows_) {
5088    useFactorization = false;
5089    factorization_->setDefaultValues();
5090
5091    int factorizationStatus = internalFactorize(0);
5092    if (factorizationStatus<0) {
5093      // some error
5094      // we should either debug or ignore
5095#ifndef NDEBUG
5096      printf("***** ClpDual strong branching factorization error - debug\n");
5097#endif
5098      return -2;
5099    } else if (factorizationStatus&&factorizationStatus<=numberRows_) {
5100      handler_->message(CLP_SINGULARITIES,messages_)
5101        <<factorizationStatus
5102        <<CoinMessageEol;
5103    }
5104  }
5105  // save stuff
5106  ClpFactorization saveFactorization(*factorization_);
5107  // save basis and solution
5108  double * saveSolution = new double[numberRows_+numberColumns_];
5109  CoinMemcpyN(solution_,
5110         numberRows_+numberColumns_,saveSolution);
5111  unsigned char * saveStatus =
5112    new unsigned char [numberRows_+numberColumns_];
5113  CoinMemcpyN(status_,numberColumns_+numberRows_,saveStatus);
5114  // save bounds as createRim makes clean copies
5115  double * saveLower = new double[numberRows_+numberColumns_];
5116  CoinMemcpyN(lower_,
5117         numberRows_+numberColumns_,saveLower);
5118  double * saveUpper = new double[numberRows_+numberColumns_];
5119  CoinMemcpyN(upper_,
5120         numberRows_+numberColumns_,saveUpper);
5121  double * saveObjective = new double[numberRows_+numberColumns_];
5122  CoinMemcpyN(cost_,
5123         numberRows_+numberColumns_,saveObjective);
5124  int * savePivot = new int [numberRows_];
5125  CoinMemcpyN(pivotVariable_, numberRows_,savePivot);
5126  // need to save/restore weights.
5127
5128  int iSolution = 0;
5129  for (i=0;i<numberVariables;i++) {
5130    int iColumn = variables[i];
5131    double objectiveChange;
5132    double saveBound;
5133
5134    // try down
5135
5136    saveBound = columnUpper_[iColumn];
5137    // external view - in case really getting optimal
5138    columnUpper_[iColumn] = newUpper[i];
5139    if (scalingFlag_<=0)
5140      upper_[iColumn] = newUpper[i]*rhsScale_;
5141    else
5142      upper_[iColumn] = (newUpper[i]/columnScale_[iColumn])*rhsScale_; // scale
5143    // Start of fast iterations
5144    int status = fastDual(alwaysFinish);
5145    CoinAssert (problemStatus_||objectiveValue_<1.0e50);
5146    // make sure plausible
5147    double obj = CoinMax(objectiveValue_,saveObjectiveValue);
5148    if (status&&problemStatus_!=3) {
5149      // not finished - might be optimal
5150      checkPrimalSolution(rowActivityWork_,columnActivityWork_);
5151      double limit = 0.0;
5152      getDblParam(ClpDualObjectiveLimit, limit);
5153      if (!numberPrimalInfeasibilities_&&obj<limit) {
5154        problemStatus_=0;
5155      }
5156      status=problemStatus_;
5157    }
5158    if (status||(problemStatus_==0&&!isDualObjectiveLimitReached())) {
5159      objectiveChange = obj-saveObjectiveValue;
5160    } else {
5161      objectiveChange = 1.0e100;
5162      status=1;
5163    }
5164    if (problemStatus_==3)
5165      status=2;
5166
5167    if (scalingFlag_<=0) {
5168      CoinMemcpyN(solution_,numberColumns_,outputSolution[iSolution]);
5169    } else {
5170      int j;
5171      double * sol = outputSolution[iSolution];
5172      for (j=0;j<numberColumns_;j++)
5173        sol[j] = solution_[j]*columnScale_[j];
5174    }
5175    outputStatus[iSolution]=status;
5176    outputIterations[iSolution]=numberIterations_;
5177    iSolution++;
5178    // restore
5179    CoinMemcpyN(saveSolution,
5180            numberRows_+numberColumns_,solution_);
5181    CoinMemcpyN(saveStatus,numberColumns_+numberRows_,status_);
5182    CoinMemcpyN(saveLower,
5183            numberRows_+numberColumns_,lower_);
5184    CoinMemcpyN(saveUpper,
5185            numberRows_+numberColumns_,upper_);
5186    CoinMemcpyN(saveObjective,
5187            numberRows_+numberColumns_,cost_);
5188    columnUpper_[iColumn] = saveBound;
5189    CoinMemcpyN(savePivot, numberRows_,pivotVariable_);
5190    delete factorization_;
5191    factorization_ = new ClpFactorization(saveFactorization);
5192
5193    newUpper[i]=objectiveChange;
5194#ifdef CLP_DEBUG
5195    printf("down on %d costs %g\n",iColumn,objectiveChange);
5196#endif
5197
5198    // try up
5199
5200    saveBound = columnLower_[iColumn];
5201    // external view - in case really getting optimal
5202    columnLower_[iColumn] = newLower[i];
5203    if (scalingFlag_<=0)
5204      lower_[iColumn] = newLower[i]*rhsScale_;
5205    else
5206      lower_[iColumn] = (newLower[i]/columnScale_[iColumn])*rhsScale_; // scale
5207    // Start of fast iterations
5208    status = fastDual(alwaysFinish);
5209    // make sure plausible
5210    obj = CoinMax(objectiveValue_,saveObjectiveValue);
5211    if (status&&problemStatus_!=3) {
5212      // not finished - might be optimal
5213      checkPrimalSolution(rowActivityWork_,columnActivityWork_);
5214      double limit = 0.0;
5215      getDblParam(ClpDualObjectiveLimit, limit);
5216      if (!numberPrimalInfeasibilities_&&obj< limit) {
5217        problemStatus_=0;
5218      }
5219      status=problemStatus_;
5220    }
5221    if (status||(problemStatus_==0&&!isDualObjectiveLimitReached())) {
5222      objectiveChange = obj-saveObjectiveValue;
5223    } else {
5224      objectiveChange = 1.0e100;
5225      status=1;
5226    }
5227    if (problemStatus_==3)
5228      status=2;
5229    if (scalingFlag_<=0) {
5230      CoinMemcpyN(solution_,numberColumns_,outputSolution[iSolution]);
5231    } else {
5232      int j;
5233      double * sol = outputSolution[iSolution];
5234      for (j=0;j<numberColumns_;j++)
5235        sol[j] = solution_[j]*columnScale_[j];
5236    }
5237    outputStatus[iSolution]=status;
5238    outputIterations[iSolution]=numberIterations_;
5239    iSolution++;
5240
5241    // restore
5242    CoinMemcpyN(saveSolution,
5243            numberRows_+numberColumns_,solution_);
5244    CoinMemcpyN(saveStatus,numberColumns_+numberRows_,status_);
5245    CoinMemcpyN(saveLower,
5246            numberRows_+numberColumns_,lower_);
5247    CoinMemcpyN(saveUpper,
5248            numberRows_+numberColumns_,upper_);
5249    CoinMemcpyN(saveObjective,
5250            numberRows_+numberColumns_,cost_);
5251    columnLower_[iColumn] = saveBound;
5252    CoinMemcpyN(savePivot, numberRows_,pivotVariable_);
5253    delete factorization_;
5254    factorization_ = new ClpFactorization(saveFactorization);
5255
5256    newLower[i]=objectiveChange;
5257#ifdef CLP_DEBUG
5258    printf("up on %d costs %g\n",iColumn,objectiveChange);
5259#endif
5260
5261    /* Possibilities are:
5262       Both sides feasible - store
5263       Neither side feasible - set objective high and exit if desired
5264       One side feasible - change bounds and resolve
5265    */
5266    if (newUpper[i]<1.0e100) {
5267      if(newLower[i]<1.0e100) {
5268        // feasible - no action
5269      } else {
5270        // up feasible, down infeasible
5271        returnCode=1;
5272        if (stopOnFirstInfeasible)
5273          break;
5274      }
5275    } else {
5276      if(newLower[i]<1.0e100) {
5277        // down feasible, up infeasible
5278        returnCode=1;
5279        if (stopOnFirstInfeasible)
5280          break;
5281      } else {
5282        // neither side feasible
5283        returnCode=-1;
5284        break;
5285      }
5286    }
5287  }
5288  delete [] saveSolution;
5289  delete [] saveLower;
5290  delete [] saveUpper;
5291  delete [] saveObjective;
5292  delete [] saveStatus;
5293  delete [] savePivot;
5294  if ((startFinishOptions&1)==0) {
5295    deleteRim(1);
5296    whatsChanged_=0;
5297  } else {
5298    // Original factorization will have been put back by last loop
5299    //delete factorization_;
5300    //factorization_ = new ClpFactorization(saveFactorization);
5301    deleteRim(0);
5302    // mark all as current
5303    whatsChanged_ = 0xffff;
5304  }
5305  objectiveValue_ = saveObjectiveValue;
5306  return returnCode;
5307}
5308// treat no pivot as finished (unless interesting)
5309int ClpSimplexDual::fastDual(bool alwaysFinish)
5310{
5311  algorithm_ = -1;
5312  secondaryStatus_=0;
5313  // Say in fast dual
5314  specialOptions_ |= 16384;
5315  int saveDont = dontFactorizePivots_;
5316  if ((specialOptions_&2048)==0)
5317    dontFactorizePivots_=0;
5318  else if(!dontFactorizePivots_)
5319    dontFactorizePivots_ = 20;
5320  //handler_->setLogLevel(63);
5321  // save data
5322  ClpDataSave data = saveData();
5323  dualTolerance_=dblParam_[ClpDualTolerance];
5324  primalTolerance_=dblParam_[ClpPrimalTolerance];
5325
5326  // save dual bound
5327  double saveDualBound = dualBound_;
5328
5329  if (alphaAccuracy_!=-1.0)
5330    alphaAccuracy_ = 1.0;
5331  double objectiveChange;
5332  // for dual we will change bounds using dualBound_
5333  // for this we need clean basis so it is after factorize
5334#if 0
5335  {
5336    int numberTotal = numberRows_+numberColumns_;
5337    double * saveSol = CoinCopyOfArray(solution_,numberTotal);
5338    double * saveDj = CoinCopyOfArray(dj_,numberTotal);
5339    double tolerance = 1.0e-8;
5340    gutsOfSolution(NULL,NULL);
5341    int j;
5342    double largestPrimal = tolerance;
5343    int iPrimal=-1;
5344    for (j=0;j<numberTotal;j++) {
5345      double difference = solution_[j]-saveSol[j];
5346      if (fabs(difference)>largestPrimal) {
5347        iPrimal=j;
5348        largestPrimal=fabs(difference);
5349      }
5350    }
5351    double largestDual = tolerance;
5352    int iDual=-1;
5353    for (j=0;j<numberTotal;j++) {
5354      double difference = dj_[j]-saveDj[j];
5355      if (fabs(difference)>largestDual&&upper_[j]>lower_[j]) {
5356        iDual=j;
5357        largestDual=fabs(difference);
5358      }
5359    }
5360    if (iPrimal>=0||iDual>=0)
5361        printf("pivots %d primal diff(%g,%d) dual diff(%g,%d)\n",
5362               factorization_->pivots(),
5363               largestPrimal,iPrimal,
5364               largestDual,iDual);
5365    delete [] saveSol;
5366    delete [] saveDj;
5367  }
5368#else
5369  if ((specialOptions_&524288)==0)
5370    gutsOfSolution(NULL,NULL);
5371#endif
5372#if 0
5373  if (numberPrimalInfeasibilities_!=1||
5374      numberDualInfeasibilities_)
5375    printf("dual %g (%d) primal %g (%d)\n",
5376           sumDualInfeasibilities_,numberDualInfeasibilities_,
5377           sumPrimalInfeasibilities_,numberPrimalInfeasibilities_);
5378#endif
5379  numberFake_ =0; // Number of variables at fake bounds
5380  numberChanged_ =0; // Number of variables with changed costs
5381  changeBounds(1,NULL,objectiveChange);
5382
5383  problemStatus_ = -1;
5384  numberIterations_=0;
5385  if ((specialOptions_&524288)==0) {
5386    factorization_->sparseThreshold(0);
5387    factorization_->goSparse();
5388  }
5389
5390  int lastCleaned=0; // last time objective or bounds cleaned up
5391
5392  // number of times we have declared optimality
5393  numberTimesOptimal_=0;
5394
5395  // This says whether to restore things etc
5396  int factorType=0;
5397  /*
5398    Status of problem:
5399    0 - optimal
5400    1 - infeasible
5401    2 - unbounded
5402    -1 - iterating
5403    -2 - factorization wanted
5404    -3 - redo checking without factorization
5405    -4 - looks infeasible
5406
5407    BUT also from whileIterating return code is:
5408
5409   -1 iterations etc
5410   -2 inaccuracy
5411   -3 slight inaccuracy (and done iterations)
5412   +0 looks optimal (might be unbounded - but we will investigate)
5413   +1 looks infeasible
5414   +3 max iterations
5415
5416  */
5417
5418  int returnCode = 0;
5419
5420  int iRow,iColumn;
5421  int maxPass = maximumIterations();
5422  while (problemStatus_<0) {
5423    // clear
5424    for (iRow=0;iRow<4;iRow++) {
5425      rowArray_[iRow]->clear();
5426    }
5427
5428    for (iColumn=0;iColumn<2;iColumn++) {
5429      columnArray_[iColumn]->clear();
5430    }
5431
5432    // give matrix (and model costs and bounds a chance to be
5433    // refreshed (normally null)
5434    matrix_->refresh(this);
5435    // If getting nowhere - why not give it a kick
5436    // does not seem to work too well - do some more work
5437    if ((specialOptions_&524288)!=0&&
5438        perturbation_<101&&numberIterations_>2*(numberRows_+numberColumns_)) {
5439      perturb();
5440      // Can't get here if values pass
5441      gutsOfSolution(NULL,NULL);
5442      if (handler_->logLevel()>2) {
5443        handler_->message(CLP_SIMPLEX_STATUS,messages_)
5444          <<numberIterations_<<objectiveValue();
5445        handler_->printing(sumPrimalInfeasibilities_>0.0)
5446          <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
5447        handler_->printing(sumDualInfeasibilities_>0.0)
5448          <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
5449        handler_->printing(numberDualInfeasibilitiesWithoutFree_
5450                           <numberDualInfeasibilities_)
5451                             <<numberDualInfeasibilitiesWithoutFree_;
5452        handler_->message()<<CoinMessageEol;
5453      }
5454    }
5455    // may factorize, checks if problem finished
5456    // should be able to speed this up on first time
5457    statusOfProblemInDual(lastCleaned,factorType,NULL,data,0);
5458
5459    // Say good factorization
5460    factorType=1;
5461    maxPass--;
5462    if (maxPass<-10) {
5463      // odd
5464      returnCode=1;
5465      problemStatus_ = 3;
5466      // can't say anything interesting - might as well return
5467#ifdef CLP_DEBUG
5468      printf("returning from fastDual after %d iterations with code %d because of loop\n",
5469             numberIterations_,returnCode);
5470#endif
5471      break;
5472    }
5473
5474    // Do iterations
5475    if (problemStatus_<0) {
5476      double * givenPi=NULL;
5477      returnCode = whileIterating(givenPi,0);
5478      if ((!alwaysFinish&&returnCode<0)||returnCode==3) {
5479        if (returnCode!=3)
5480          assert (problemStatus_<0);
5481        returnCode=1;
5482        problemStatus_ = 3;
5483        // can't say anything interesting - might as well return
5484#ifdef CLP_DEBUG
5485        printf("returning from fastDual after %d iterations with code %d\n",
5486               numberIterations_,returnCode);
5487#endif
5488        break;
5489      }
5490      if (returnCode==-2)
5491        factorType=3;
5492      returnCode=0;
5493    }
5494  }
5495
5496  // clear
5497  for (iRow=0;iRow<4;iRow++) {
5498    rowArray_[iRow]->clear();
5499  }
5500
5501  for (iColumn=0;iColumn<2;iColumn++) {
5502    columnArray_[iColumn]->clear();
5503  }
5504  // Say not in fast dual
5505  specialOptions_ &= ~16384;
5506  assert(!numberFake_||((specialOptions_&(2048|4096))!=0&&dualBound_>=1.0e8)
5507         ||returnCode||problemStatus_); // all bounds should be okay
5508  if (numberFake_>0) {
5509    // Set back
5510    double dummy;
5511    changeBounds(2,NULL,dummy);
5512  }
5513  // Restore any saved stuff
5514  restoreData(data);
5515  dontFactorizePivots_ = saveDont;
5516  dualBound_ = saveDualBound;
5517  return returnCode;
5518}
5519// This does first part of StrongBranching
5520ClpFactorization *
5521ClpSimplexDual::setupForStrongBranching(char * arrays, int numberRows, int numberColumns)
5522{
5523  algorithm_ = -1;
5524  // put in standard form (and make row copy)
5525  // create modifiable copies of model rim and do optional scaling
5526  int startFinishOptions;
5527  /*  COIN_CLP_VETTED
5528      Looks safe for Cbc
5529  */
5530  bool safeOption = ((specialOptions_&COIN_CBC_USING_CLP)!=0);
5531  if((specialOptions_&4096)==0||(auxiliaryModel_&&!safeOption)) {
5532    startFinishOptions=0;
5533  } else {
5534    startFinishOptions=1+2+4;
5535  }
5536  createRim(7+8+16+32,true,startFinishOptions);
5537  // Do initial factorization
5538  // and set certain stuff
5539  // We can either set increasing rows so ...IsBasic gives pivot row
5540  // or we can just increment iBasic one by one
5541  // for now let ...iBasic give pivot row
5542  bool useFactorization=false;
5543  if ((startFinishOptions&2)!=0&&(whatsChanged_&(2+512))==2+512) {
5544    useFactorization=true; // Keep factorization if possible
5545    // switch off factorization if bad
5546    if (pivotVariable_[0]<0||factorization_->numberRows()!=numberRows_)
5547      useFactorization=false;
5548  }
5549  if (!useFactorization) {
5550    factorization_->setDefaultValues();
5551
5552    int factorizationStatus = internalFactorize(0);
5553    if (factorizationStatus<0) {
5554      // some error
5555      // we should either debug or ignore
5556#ifndef NDEBUG
5557      printf("***** ClpDual strong branching factorization error - debug\n");
5558#endif
5559    } else if (factorizationStatus&&factorizationStatus<=numberRows_) {
5560      handler_->message(CLP_SINGULARITIES,messages_)
5561        <<factorizationStatus
5562        <<CoinMessageEol;
5563    }
5564  }
5565  double * arrayD = (double *) arrays;
5566  arrayD[0]=objectiveValue()*optimizationDirection_;
5567  double * saveSolution = arrayD+1;
5568  double * saveLower = saveSolution + (numberRows+numberColumns);
5569  double * saveUpper = saveLower + (numberRows+numberColumns);
5570  double * saveObjective = saveUpper + (numberRows+numberColumns);
5571  double * saveLowerOriginal = saveObjective + (numberRows+numberColumns);
5572  double * saveUpperOriginal = saveLowerOriginal + numberColumns;
5573  arrayD = saveUpperOriginal + numberColumns;
5574  int * savePivot = (int *) arrayD;
5575  int * whichRow = savePivot+numberRows;
5576  int * whichColumn = whichRow + 3*numberRows;
5577  int * arrayI = whichColumn + 2*numberColumns;
5578  unsigned char * saveStatus = (unsigned char *) (arrayI+1);
5579  // save stuff
5580  // save basis and solution
5581  CoinMemcpyN(solution_,
5582          numberRows_+numberColumns_,saveSolution);
5583  CoinMemcpyN(status_,numberColumns_+numberRows_,saveStatus);
5584  CoinMemcpyN(lower_,
5585          numberRows_+numberColumns_,saveLower);
5586  CoinMemcpyN(upper_,
5587          numberRows_+numberColumns_,saveUpper);
5588  CoinMemcpyN(cost_,
5589          numberRows_+numberColumns_,saveObjective);
5590  CoinMemcpyN(pivotVariable_, numberRows_,savePivot);
5591  return new ClpFactorization(*factorization_);
5592}
5593// This cleans up after strong branching
5594void
5595ClpSimplexDual::cleanupAfterStrongBranching()
5596{
5597  int startFinishOptions;
5598  /*  COIN_CLP_VETTED
5599      Looks safe for Cbc
5600  */
5601  bool safeOption = ((specialOptions_&COIN_CBC_USING_CLP)!=0);
5602  if((specialOptions_&4096)==0||(auxiliaryModel_&&!safeOption)) {
5603    startFinishOptions=0;
5604  } else {
5605    startFinishOptions=1+2+4;
5606  }
5607  if ((startFinishOptions&1)==0) {
5608    deleteRim(1);
5609    whatsChanged_=0;
5610  } else {
5611    // Original factorization will have been put back by last loop
5612    //delete factorization_;
5613    //factorization_ = new ClpFactorization(saveFactorization);
5614    deleteRim(0);
5615    // mark all as current
5616    whatsChanged_ = 0xffff;
5617  }
5618}
5619/* Checks number of variables at fake bounds.  This is used by fastDual
5620   so can exit gracefully before end */
5621int
5622ClpSimplexDual::numberAtFakeBound()
5623{
5624  int iSequence;
5625  int numberFake=0;
5626
5627  for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
5628    FakeBound bound = getFakeBound(iSequence);
5629    switch(getStatus(iSequence)) {
5630
5631    case basic:
5632      break;
5633    case isFree:
5634    case superBasic:
5635    case ClpSimplex::isFixed:
5636      //setFakeBound (iSequence, noFake);
5637      break;
5638    case atUpperBound:
5639      if (bound==upperFake||bound==bothFake)
5640        numberFake++;
5641      break;
5642    case atLowerBound:
5643      if (bound==lowerFake||bound==bothFake)
5644        numberFake++;
5645      break;
5646    }
5647  }
5648  numberFake_ = numberFake;
5649  return numberFake;
5650}
5651/* Pivot out a variable and choose an incoing one.  Assumes dual
5652   feasible - will not go through a reduced cost.
5653   Returns step length in theta
5654   Returns ray in ray_ (or NULL if no pivot)
5655   Return codes as before but -1 means no acceptable pivot
5656*/
5657int
5658ClpSimplexDual::pivotResult()
5659{
5660  abort();
5661  return 0;
5662}
5663/*
5664   Row array has row part of pivot row
5665   Column array has column part.
5666   This is used in dual values pass
5667*/
5668void
5669ClpSimplexDual::checkPossibleValuesMove(CoinIndexedVector * rowArray,
5670                                        CoinIndexedVector * columnArray,
5671                                        double acceptablePivot)
5672{
5673  double * work;
5674  int number;
5675  int * which;
5676  int iSection;
5677
5678  double tolerance = dualTolerance_*1.001;
5679
5681  double changeDown ;
5682  double thetaUp = 1.0e31;
5684  double bestAlphaUp = acceptablePivot*0.99999;
5685  int sequenceDown =-1;
5686  int sequenceUp = sequenceOut_;
5687
5688  double djBasic = dj_[sequenceOut_];
5689  if (djBasic>0.0) {
5690    // basic at lower bound so directionOut_ 1 and -1 in pivot row
5691    // dj will go to zero on other way
5692    thetaUp = djBasic;
5693    changeDown = -lower_[sequenceOut_];
5694  } else {
5695    // basic at upper bound so directionOut_ -1 and 1 in pivot row
5696    // dj will go to zero on other way
5697    thetaUp = -djBasic;
5698    changeDown = upper_[sequenceOut_];
5699  }
5700  bestAlphaUp = 1.0;
5702
5703  double alphaUp=0.0;
5705
5706  for (iSection=0;iSection<2;iSection++) {
5707
5708    int i;
5709    if (!iSection) {
5710      work = rowArray->denseVector();
5711      number = rowArray->getNumElements();
5712      which = rowArray->getIndices();
5714    } else {
5715      work = columnArray->denseVector();
5716      number = columnArray->getNumElements();
5717      which = columnArray->getIndices();
5719    }
5720
5721    for (i=0;i<number;i++) {
5722      int iSequence = which[i];
5723      int iSequence2 = iSequence + addSequence;
5724      double alpha;
5725      double oldValue;
5726      double value;
5727
5728      switch(getStatus(iSequence2)) {
5729
5730      case basic:
5731        break;
5732      case ClpSimplex::isFixed:
5733        alpha = work[i];
5734        changeDown += alpha*upper_[iSequence2];
5735        break;
5736      case isFree:
5737      case superBasic:
5738        alpha = work[i];
5739        // dj must be effectively zero as dual feasible
5740        if (fabs(alpha)>bestAlphaUp) {
5742          thetaUp = 0.0;
5745          sequenceDown =iSequence2;
5746          sequenceUp = sequenceDown;
5747          alphaUp = alpha;
5749        }
5750        break;
5751      case atUpperBound:
5752        alpha = work[i];
5753        oldValue = dj_[iSequence2];
5754        changeDown += alpha*upper_[iSequence2];
5755        if (alpha>=acceptablePivot) {
5756          // might do other way
5757          value = oldValue+thetaUp*alpha;
5758          if (value>-tolerance) {
5759            if (value>tolerance||fabs(alpha)>bestAlphaUp) {
5760              thetaUp = -oldValue/alpha;
5761              bestAlphaUp = fabs(alpha);
5762              sequenceUp = iSequence2;
5763              alphaUp = alpha;
5764            }
5765          }
5766        } else if (alpha<=-acceptablePivot) {
5767          // might do this way
5769          if (value>-tolerance) {
5773              sequenceDown = iSequence2;
5775            }
5776          }
5777        }
5778        break;
5779      case atLowerBound:
5780        alpha = work[i];
5781        oldValue = dj_[iSequence2];
5782        changeDown += alpha*lower_[iSequence2];
5783        if (alpha<=-acceptablePivot) {
5784          // might do other way
5785          value = oldValue+thetaUp*alpha;
5786          if (value<tolerance) {
5787            if (value<-tolerance||fabs(alpha)>bestAlphaUp) {
5788              thetaUp = -oldValue/alpha;
5789              bestAlphaUp = fabs(alpha);
5790              sequenceUp = iSequence2;
5791              alphaUp = alpha;
5792            }
5793          }
5794        } else if (alpha>=acceptablePivot) {
5795          // might do this way
5797          if (value<tolerance) {
5801              sequenceDown = iSequence2;
5803            }
5804          }
5805        }
5806        break;
5807      }
5808    }
5809  }
5810  thetaUp *= -1.0;
5811  double changeUp = -thetaUp * changeDown;
5814    // largest
5816      sequenceDown =-1;
5817    }
5818  }
5819  // choose
5820  sequenceIn_=-1;
5821  if (changeDown>changeUp&&sequenceDown>=0) {
5823    if (fabs(changeDown)<1.0e30)
5824      sequenceIn_ = sequenceDown;
5826#ifdef CLP_DEBUG
5827    if ((handler_->logLevel()&32))
5828      printf("predicted way - dirout %d, change %g,%g theta %g\n",
5829             directionOut_,changeDown,changeUp,theta_);
5830#endif
5831  } else {
5832    theta_ = thetaUp;
5833    if (fabs(changeUp)<1.0e30)
5834      sequenceIn_ = sequenceUp;
5835    alpha_ = alphaUp;
5836    if (sequenceIn_!=sequenceOut_) {
5837#ifdef CLP_DEBUG
5838      if ((handler_->logLevel()&32))
5839        printf("opposite way - dirout %d, change %g,%g theta %g\n",
5840               directionOut_,changeDown,changeUp,theta_);
5841#endif
5842    } else {
5843#ifdef CLP_DEBUG
5844      if ((handler_->logLevel()&32))
5845        printf("opposite way to zero dj - dirout %d, change %g,%g theta %g\n",
5846               directionOut_,changeDown,changeUp,theta_);
5847#endif
5848    }
5849  }
5850  if (sequenceIn_>=0) {
5851    lowerIn_ = lower_[sequenceIn_];
5852    upperIn_ = upper_[sequenceIn_];
5853    valueIn_ = solution_[sequenceIn_];
5854    dualIn_ = dj_[sequenceIn_];
5855
5856    if (alpha_<0.0) {
5857      // as if from upper bound
5858      directionIn_=-1;
5859      upperIn_=valueIn_;
5860    } else {
5861      // as if from lower bound
5862      directionIn_=1;
5863      lowerIn_=valueIn_;
5864    }
5865  }
5866}
5867/*
5868   Row array has row part of pivot row
5869   Column array has column part.
5870   This is used in cleanup
5871*/
5872void
5873ClpSimplexDual::checkPossibleCleanup(CoinIndexedVector * rowArray,
5874                                        CoinIndexedVector * columnArray,
5875                                        double acceptablePivot)
5876{
5877  double * work;
5878  int number;
5879  int * which;
5880  int iSection;
5881
5882  double tolerance = dualTolerance_*1.001;
5883