source: trunk/Clp/src/ClpSimplexDual.cpp @ 2063

Last change on this file since 2063 was 2063, checked in by forrest, 5 years ago

fix typo and add some stability fixes

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