# source:stable/1.10/Clp/src/ClpSimplexDual.cpp@1404

Last change on this file since 1404 was 1404, checked in by tkr, 11 years ago

Merging r1375:1393 from stable/BSP to stable/1.10

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