# source:trunk/ClpSimplexDual.cpp@659

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

fix bug

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