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

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

improve handling of free variables in dual

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