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

Last change on this file since 1723 was 1723, checked in by forrest, 9 years ago

out some printf statements

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