source: trunk/ClpSimplexDual.cpp @ 653

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

for slim

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