source: stable/1.16/Clp/src/ClpSimplexDual.cpp @ 2114

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

was not totally threadsafe!

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 337.8 KB
Line 
1/* $Id: ClpSimplexDual.cpp 2114 2015-02-10 12:12:46Z forrest $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6
7/* Notes on implementation of dual simplex algorithm.
8
9   When dual feasible:
10
11   If primal feasible, we are optimal.  Otherwise choose an infeasible
12   basic variable to leave basis (normally going to nearest bound) (B).  We
13   now need to find an incoming variable which will leave problem
14   dual feasible so we get the row of the tableau corresponding to
15   the basic variable (with the correct sign depending if basic variable
16   above or below feasibility region - as that affects whether reduced
17   cost on outgoing variable has to be positive or negative).
18
19   We now perform a ratio test to determine which incoming variable will
20   preserve dual feasibility (C).  If no variable found then problem
21   is infeasible (in primal sense).  If there is a variable, we then
22   perform pivot and repeat.  Trivial?
23
24   -------------------------------------------
25
26   A) How do we get dual feasible?  If all variables have bounds then
27   it is trivial to get feasible by putting non-basic variables to
28   correct bounds.  OSL did not have a phase 1/phase 2 approach but
29   instead effectively put fake bounds on variables and this is the
30   approach here, although I had hoped to make it cleaner.
31
32   If there is a weight of X on getting dual feasible:
33     Non-basic variables with negative reduced costs are put to
34     lesser of their upper bound and their lower bound + X.
35     Similarly, mutatis mutandis, for positive reduced costs.
36
37   Free variables should normally be in basis, otherwise I have
38   coding which may be able to come out (and may not be correct).
39
40   In OSL, this weight was changed heuristically, here at present
41   it is only increased if problem looks finished.  If problem is
42   feasible I check for unboundedness.  If not unbounded we
43   could play with going into primal.  As long as weights increase
44   any algorithm would be finite.
45
46   B) Which outgoing variable to choose is a virtual base class.
47   For difficult problems steepest edge is preferred while for
48   very easy (large) problems we will need partial scan.
49
50   C) Sounds easy, but this is hardest part of algorithm.
51      1) Instead of stopping at first choice, we may be able
52      to flip that variable to other bound and if objective
53      still improving choose again.  These mini iterations can
54      increase speed by orders of magnitude but we may need to
55      go to more of a bucket choice of variable rather than looking
56      at them one by one (for speed).
57      2) Accuracy.  Reduced costs may be of wrong sign but less than
58      tolerance.  Pivoting on these makes objective go backwards.
59      OSL modified cost so a zero move was made, Gill et al
60      (in primal analogue) modified so a strictly positive move was
61      made.  It is not quite as neat in dual but that is what we
62      try and do.  The two problems are that re-factorizations can
63      change reduced costs above and below tolerances and that when
64      finished we need to reset costs and try again.
65      3) Degeneracy.  Gill et al helps but may not be enough.  We
66      may need more.  Also it can improve speed a lot if we perturb
67      the costs significantly.
68
69  References:
70     Forrest and Goldfarb, Steepest-edge simplex algorithms for
71       linear programming - Mathematical Programming 1992
72     Forrest and Tomlin, Implementing the simplex method for
73       the Optimization Subroutine Library - IBM Systems Journal 1992
74     Gill, Murray, Saunders, Wright A Practical Anti-Cycling
75       Procedure for Linear and Nonlinear Programming SOL report 1988
76
77
78  TODO:
79
80  a) Better recovery procedures.  At present I never check on forward
81     progress.  There is checkpoint/restart with reducing
82     re-factorization frequency, but this is only on singular
83     factorizations.
84  b) Fast methods for large easy problems (and also the option for
85     the code to automatically choose which method).
86  c) We need to be able to stop in various ways for OSI - this
87     is fairly easy.
88
89 */
90#ifdef COIN_DEVELOP
91#undef COIN_DEVELOP
92#define COIN_DEVELOP 2
93#endif
94
95#include "CoinPragma.hpp"
96
97#include <math.h>
98
99#include "CoinHelperFunctions.hpp"
100#include "ClpHelperFunctions.hpp"
101#include "ClpSimplexDual.hpp"
102#include "ClpEventHandler.hpp"
103#include "ClpFactorization.hpp"
104#include "CoinPackedMatrix.hpp"
105#include "CoinIndexedVector.hpp"
106#include "CoinFloatEqual.hpp"
107#include "ClpDualRowDantzig.hpp"
108#include "ClpMessage.hpp"
109#include "ClpLinearObjective.hpp"
110#include <cfloat>
111#include <cassert>
112#include <string>
113#include <stdio.h>
114#include <iostream>
115//#define CLP_DEBUG 1
116// To force to follow another run put logfile name here and define
117//#define FORCE_FOLLOW
118#ifdef FORCE_FOLLOW
119static FILE * fpFollow = NULL;
120static char * forceFile = "old.log";
121static int force_in = -1;
122static int force_out = -1;
123static int force_iteration = 0;
124#endif
125//#define VUB
126#ifdef VUB
127extern int * vub;
128extern int * toVub;
129extern int * nextDescendent;
130#endif
131#ifdef NDEBUG
132#define NDEBUG_CLP
133#endif
134#ifndef CLP_INVESTIGATE
135#define NDEBUG_CLP
136#endif
137// dual
138
139/* *** Method
140   This is a vanilla version of dual simplex.
141
142   It tries to be a single phase approach with a weight of 1.0 being
143   given to getting optimal and a weight of dualBound_ being
144   given to getting dual feasible.  In this version I have used the
145   idea that this weight can be thought of as a fake bound.  If the
146   distance between the lower and upper bounds on a variable is less
147   than the feasibility weight then we are always better off flipping
148   to other bound to make dual feasible.  If the distance is greater
149   then we make up a fake bound dualBound_ away from one bound.
150   If we end up optimal or primal infeasible, we check to see if
151   bounds okay.  If so we have finished, if not we increase dualBound_
152   and continue (after checking if unbounded). I am undecided about
153   free variables - there is coding but I am not sure about it.  At
154   present I put them in basis anyway.
155
156   The code is designed to take advantage of sparsity so arrays are
157   seldom zeroed out from scratch or gone over in their entirety.
158   The only exception is a full scan to find outgoing variable.  This
159   will be changed to keep an updated list of infeasibilities (or squares
160   if steepest edge).  Also on easy problems we don't need full scan - just
161   pick first reasonable.
162
163   One problem is how to tackle degeneracy and accuracy.  At present
164   I am using the modification of costs which I put in OSL and which was
165   extended by Gill et al.  I am still not sure of the exact details.
166
167   The flow of dual is three while loops as follows:
168
169   while (not finished) {
170
171     while (not clean solution) {
172
173        Factorize and/or clean up solution by flipping variables so
174  dual feasible.  If looks finished check fake dual bounds.
175  Repeat until status is iterating (-1) or finished (0,1,2)
176
177     }
178
179     while (status==-1) {
180
181       Iterate until no pivot in or out or time to re-factorize.
182
183       Flow is:
184
185       choose pivot row (outgoing variable).  if none then
186 we are primal feasible so looks as if done but we need to
187 break and check bounds etc.
188
189 Get pivot row in tableau
190
191       Choose incoming column.  If we don't find one then we look
192 primal infeasible so break and check bounds etc.  (Also the
193 pivot tolerance is larger after any iterations so that may be
194 reason)
195
196       If we do find incoming column, we may have to adjust costs to
197 keep going forwards (anti-degeneracy).  Check pivot will be stable
198 and if unstable throw away iteration (we will need to implement
199 flagging of basic variables sometime) and break to re-factorize.
200 If minor error re-factorize after iteration.
201
202 Update everything (this may involve flipping variables to stay
203 dual feasible.
204
205     }
206
207   }
208
209   At present we never check we are going forwards.  I overdid that in
210   OSL so will try and make a last resort.
211
212   Needs partial scan pivot out option.
213   Needs dantzig, uninitialized and full steepest edge options (can still
214   use partial scan)
215
216   May need other anti-degeneracy measures, especially if we try and use
217   loose tolerances as a way to solve in fewer iterations.
218
219   I like idea of dynamic scaling.  This gives opportunity to decouple
220   different implications of scaling for accuracy, iteration count and
221   feasibility tolerance.
222
223*/
224#define CLEAN_FIXED 0
225// Startup part of dual (may be extended to other algorithms)
226int
227ClpSimplexDual::startupSolve(int ifValuesPass, double * saveDuals, int startFinishOptions)
228{
229     // If values pass then save given duals round check solution
230     // sanity check
231     // initialize - no values pass and algorithm_ is -1
232     // put in standard form (and make row copy)
233     // create modifiable copies of model rim and do optional scaling
234     // If problem looks okay
235     // Do initial factorization
236     // If user asked for perturbation - do it
237     numberFake_ = 0; // Number of variables at fake bounds
238     numberChanged_ = 0; // Number of variables with changed costs
239     if (!startup(0, startFinishOptions)) {
240          int usePrimal = 0;
241          // looks okay
242          // Superbasic variables not allowed
243          // If values pass then scale pi
244          if (ifValuesPass) {
245               if (problemStatus_ && perturbation_ < 100)
246                    usePrimal = perturb();
247               int i;
248               if (scalingFlag_ > 0) {
249                    for (i = 0; i < numberRows_; i++) {
250                         dual_[i] = saveDuals[i] * inverseRowScale_[i];
251                    }
252               } else {
253                    CoinMemcpyN(saveDuals, numberRows_, dual_);
254               }
255               // now create my duals
256               for (i = 0; i < numberRows_; i++) {
257                    // slack
258                    double value = dual_[i];
259                    value += rowObjectiveWork_[i];
260                    saveDuals[i+numberColumns_] = value;
261               }
262               CoinMemcpyN(objectiveWork_, numberColumns_, saveDuals);
263               transposeTimes(-1.0, dual_, saveDuals);
264               // make reduced costs okay
265               for (i = 0; i < numberColumns_; i++) {
266                    if (getStatus(i) == atLowerBound) {
267                         if (saveDuals[i] < 0.0) {
268                              //if (saveDuals[i]<-1.0e-3)
269                              //printf("bad dj at lb %d %g\n",i,saveDuals[i]);
270                              saveDuals[i] = 0.0;
271                         }
272                    } else if (getStatus(i) == atUpperBound) {
273                         if (saveDuals[i] > 0.0) {
274                              //if (saveDuals[i]>1.0e-3)
275                              //printf("bad dj at ub %d %g\n",i,saveDuals[i]);
276                              saveDuals[i] = 0.0;
277                         }
278                    }
279               }
280               CoinMemcpyN(saveDuals, (numberColumns_ + numberRows_), dj_);
281               // set up possible ones
282               for (i = 0; i < numberRows_ + numberColumns_; i++)
283                    clearPivoted(i);
284               int iRow;
285               for (iRow = 0; iRow < numberRows_; iRow++) {
286                    int iPivot = pivotVariable_[iRow];
287                    if (fabs(saveDuals[iPivot]) > dualTolerance_) {
288                         if (getStatus(iPivot) != isFree)
289                              setPivoted(iPivot);
290                    }
291               }
292          } else if ((specialOptions_ & 1024) != 0 && CLEAN_FIXED) {
293               // set up possible ones
294               for (int i = 0; i < numberRows_ + numberColumns_; i++)
295                    clearPivoted(i);
296               int iRow;
297               for (iRow = 0; iRow < numberRows_; iRow++) {
298                    int iPivot = pivotVariable_[iRow];
299                    if (iPivot < numberColumns_ && lower_[iPivot] == upper_[iPivot]) {
300                         setPivoted(iPivot);
301                    }
302               }
303          }
304
305          double objectiveChange;
306          assert (!numberFake_);
307          assert (numberChanged_ == 0);
308          if (!numberFake_) // if nonzero then adjust
309               changeBounds(1, NULL, objectiveChange);
310
311          if (!ifValuesPass) {
312               // Check optimal
313               if (!numberDualInfeasibilities_ && !numberPrimalInfeasibilities_)
314                    problemStatus_ = 0;
315          }
316          if (problemStatus_ < 0 && perturbation_ < 100) {
317               bool inCbcOrOther = (specialOptions_ & 0x03000000) != 0;
318               if (!inCbcOrOther)
319                    usePrimal = perturb();
320               // Can't get here if values pass
321               gutsOfSolution(NULL, NULL);
322#ifdef CLP_INVESTIGATE
323               if (numberDualInfeasibilities_)
324                    printf("ZZZ %d primal %d dual - sumdinf %g\n",
325                           numberPrimalInfeasibilities_,
326                           numberDualInfeasibilities_, sumDualInfeasibilities_);
327#endif
328               if (handler_->logLevel() > 2) {
329                    handler_->message(CLP_SIMPLEX_STATUS, messages_)
330                              << numberIterations_ << objectiveValue();
331                    handler_->printing(sumPrimalInfeasibilities_ > 0.0)
332                              << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
333                    handler_->printing(sumDualInfeasibilities_ > 0.0)
334                              << sumDualInfeasibilities_ << numberDualInfeasibilities_;
335                    handler_->printing(numberDualInfeasibilitiesWithoutFree_
336                                       < numberDualInfeasibilities_)
337                              << numberDualInfeasibilitiesWithoutFree_;
338                    handler_->message() << CoinMessageEol;
339               }
340               if (inCbcOrOther) {
341                    if (numberPrimalInfeasibilities_) {
342                         usePrimal = perturb();
343                         if (perturbation_ >= 101) {
344                              computeDuals(NULL);
345                              //gutsOfSolution(NULL,NULL);
346                              checkDualSolution(); // recompute objective
347                         }
348                    } else if (numberDualInfeasibilities_) {
349                         problemStatus_ = 10;
350                         if ((moreSpecialOptions_ & 32) != 0 && false)
351                              problemStatus_ = 0; // say optimal!!
352#if COIN_DEVELOP>2
353
354                         printf("returning at %d\n", __LINE__);
355#endif
356                         return 1; // to primal
357                    }
358               }
359          } else if (!ifValuesPass) {
360               gutsOfSolution(NULL, NULL);
361               // double check
362               if (numberDualInfeasibilities_ || numberPrimalInfeasibilities_)
363                    problemStatus_ = -1;
364          }
365          if (usePrimal) {
366               problemStatus_ = 10;
367#if COIN_DEVELOP>2
368               printf("returning to use primal (no obj) at %d\n", __LINE__);
369#endif
370          }
371          return usePrimal;
372     } else {
373          return 1;
374     }
375}
376void
377ClpSimplexDual::finishSolve(int startFinishOptions)
378{
379     assert (problemStatus_ || !sumPrimalInfeasibilities_);
380
381     // clean up
382     finish(startFinishOptions);
383}
384//#define CLP_REPORT_PROGRESS
385#ifdef CLP_REPORT_PROGRESS
386static int ixxxxxx = 0;
387static int ixxyyyy = 90;
388#endif
389#ifdef CLP_INVESTIGATE_SERIAL
390static int z_reason[7] = {0, 0, 0, 0, 0, 0, 0};
391static int z_thinks = -1;
392#endif
393void
394ClpSimplexDual::gutsOfDual(int ifValuesPass, double * & saveDuals, int initialStatus,
395                           ClpDataSave & data)
396{
397#ifdef CLP_INVESTIGATE_SERIAL
398     z_reason[0]++;
399     z_thinks = -1;
400     int nPivots = 9999;
401#endif
402     double largestPrimalError = 0.0;
403     double largestDualError = 0.0;
404     int numberRayTries=0;
405     // Start can skip some things in transposeTimes
406     specialOptions_ |= 131072;
407     int lastCleaned = 0; // last time objective or bounds cleaned up
408
409     // This says whether to restore things etc
410     // startup will have factorized so can skip
411     int factorType = 0;
412     // Start check for cycles
413     progress_.startCheck();
414     // Say change made on first iteration
415     changeMade_ = 1;
416     // Say last objective infinite
417     //lastObjectiveValue_=-COIN_DBL_MAX;
418     progressFlag_ = 0;
419     /*
420       Status of problem:
421       0 - optimal
422       1 - infeasible
423       2 - unbounded
424       -1 - iterating
425       -2 - factorization wanted
426       -3 - redo checking without factorization
427       -4 - looks infeasible
428     */
429     while (problemStatus_ < 0) {
430          int iRow, iColumn;
431          // clear
432          for (iRow = 0; iRow < 4; iRow++) {
433               rowArray_[iRow]->clear();
434          }
435
436          for (iColumn = 0; iColumn < 2; iColumn++) {
437               columnArray_[iColumn]->clear();
438          }
439
440          // give matrix (and model costs and bounds a chance to be
441          // refreshed (normally null)
442          matrix_->refresh(this);
443          // If getting nowhere - why not give it a kick
444          // does not seem to work too well - do some more work
445          if (perturbation_ < 101 && numberIterations_ > 2 * (numberRows_ + numberColumns_) && (moreSpecialOptions_&1048576)==0
446                    && initialStatus != 10) {
447               perturb();
448               // Can't get here if values pass
449               gutsOfSolution(NULL, NULL);
450               if (handler_->logLevel() > 2) {
451                    handler_->message(CLP_SIMPLEX_STATUS, messages_)
452                              << numberIterations_ << objectiveValue();
453                    handler_->printing(sumPrimalInfeasibilities_ > 0.0)
454                              << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
455                    handler_->printing(sumDualInfeasibilities_ > 0.0)
456                              << sumDualInfeasibilities_ << numberDualInfeasibilities_;
457                    handler_->printing(numberDualInfeasibilitiesWithoutFree_
458                                       < numberDualInfeasibilities_)
459                              << numberDualInfeasibilitiesWithoutFree_;
460                    handler_->message() << CoinMessageEol;
461               }
462          }
463          // see if in Cbc etc
464          bool inCbcOrOther = (specialOptions_ & 0x03000000) != 0;
465#if 0
466          bool gotoPrimal = false;
467          if (inCbcOrOther && numberIterations_ > disasterArea_ + numberRows_ &&
468                    numberDualInfeasibilitiesWithoutFree_ && largestDualError_ > 1.0e-1) {
469               if (!disasterArea_) {
470                    printf("trying all slack\n");
471                    // try all slack basis
472                    allSlackBasis(true);
473                    disasterArea_ = 2 * numberRows_;
474               } else {
475                    printf("going to primal\n");
476                    // go to primal
477                    gotoPrimal = true;
478                    allSlackBasis(true);
479               }
480          }
481#endif
482          bool disaster = false;
483          if (disasterArea_ && inCbcOrOther && disasterArea_->check()) {
484               disasterArea_->saveInfo();
485               disaster = true;
486          }
487          // may factorize, checks if problem finished
488          statusOfProblemInDual(lastCleaned, factorType, saveDuals, data,
489                                ifValuesPass);
490          if ((specialOptions_&2097152)!=0&&problemStatus_==1&&!ray_&&
491              !numberRayTries) {
492            numberRayTries=1;
493            problemStatus_=-1;
494          }
495          largestPrimalError = CoinMax(largestPrimalError, largestPrimalError_);
496          largestDualError = CoinMax(largestDualError, largestDualError_);
497          if (disaster)
498               problemStatus_ = 3;
499          // If values pass then do easy ones on first time
500          if (ifValuesPass &&
501                    progress_.lastIterationNumber(0) < 0 && saveDuals) {
502               doEasyOnesInValuesPass(saveDuals);
503          }
504
505          // Say good factorization
506          factorType = 1;
507          if (data.sparseThreshold_) {
508               // use default at present
509               factorization_->sparseThreshold(0);
510               factorization_->goSparse();
511          }
512
513          // exit if victory declared
514          if (problemStatus_ >= 0)
515               break;
516
517          // test for maximum iterations
518          if (hitMaximumIterations() || (ifValuesPass == 2 && !saveDuals)) {
519               problemStatus_ = 3;
520               break;
521          }
522          if (ifValuesPass && !saveDuals) {
523               // end of values pass
524               ifValuesPass = 0;
525               int status = eventHandler_->event(ClpEventHandler::endOfValuesPass);
526               if (status >= 0) {
527                    problemStatus_ = 5;
528                    secondaryStatus_ = ClpEventHandler::endOfValuesPass;
529                    break;
530               }
531          }
532          // Check event
533          {
534               int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
535               if (status >= 0) {
536                    problemStatus_ = 5;
537                    secondaryStatus_ = ClpEventHandler::endOfFactorization;
538                    break;
539               }
540          }
541          // 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                    if (value < -tolerance) {
2336                         // flipping bounds
2337                         double movement = mult * (lower[iSequence] - upper[iSequence]);
2338                         which[numberInfeasibilities++] = iSequence;
2339#ifndef NDEBUG
2340                         if (fabs(movement) >= 1.0e30)
2341                              resetFakeBounds(-1000 - iSequence);
2342#endif
2343#ifdef CLP_DEBUG
2344                         if ((handler_->logLevel() & 32))
2345                              printf("%d %d, new dj %g, alpha %g, movement %g\n",
2346                                     0, iSequence, value, alphaI, movement);
2347#endif
2348                         changeObj -= movement * cost[iSequence];
2349                         outputArray->quickAdd(iSequence, movement);
2350                    }
2351               }
2352          }
2353          // Do columns
2354          reducedCost = djRegion(1);
2355          lower = lowerRegion(1);
2356          upper = upperRegion(1);
2357          cost = costRegion(1);
2358          // set number of infeasibilities in row array
2359          numberRowInfeasibilities = numberInfeasibilities;
2360          rowArray->setNumElements(numberInfeasibilities);
2361          numberInfeasibilities = 0;
2362          work = columnArray->denseVector();
2363          number = columnArray->getNumElements();
2364          which = columnArray->getIndices();
2365          if ((moreSpecialOptions_ & 8) != 0) {
2366               const unsigned char * COIN_RESTRICT statusArray = status_;
2367               for (i = 0; i < number; i++) {
2368                    int iSequence = which[i];
2369                    double alphaI = work[i];
2370                    work[i] = 0.0;
2371
2372                    int iStatus = (statusArray[iSequence] & 3) - 1;
2373                    if (iStatus) {
2374                         double value = reducedCost[iSequence] - theta * alphaI;
2375                         reducedCost[iSequence] = value;
2376                         double mult = multiplier[iStatus-1];
2377                         value *= mult;
2378                         if (value < -tolerance) {
2379                              // flipping bounds
2380                              double movement = mult * (upper[iSequence] - lower[iSequence]);
2381                              which[numberInfeasibilities++] = iSequence;
2382#ifndef NDEBUG
2383                              if (fabs(movement) >= 1.0e30)
2384                                   resetFakeBounds(-1000 - iSequence);
2385#endif
2386#ifdef CLP_DEBUG
2387                              if ((handler_->logLevel() & 32))
2388                                   printf("%d %d, new dj %g, alpha %g, movement %g\n",
2389                                          1, iSequence, value, alphaI, movement);
2390#endif
2391                              changeObj += movement * cost[iSequence];
2392                              matrix_->add(this, outputArray, iSequence, movement);
2393                         }
2394                    }
2395               }
2396          } else {
2397               for (i = 0; i < number; i++) {
2398                    int iSequence = which[i];
2399                    double alphaI = work[i];
2400                    work[i] = 0.0;
2401
2402                    Status status = getStatus(iSequence);
2403                    if (status == atLowerBound) {
2404                         double value = reducedCost[iSequence] - theta * alphaI;
2405                         reducedCost[iSequence] = value;
2406                         double movement = 0.0;
2407
2408                         if (value < -tolerance) {
2409                              // to upper bound
2410                              which[numberInfeasibilities++] = iSequence;
2411                              movement = upper[iSequence] - lower[iSequence];
2412#ifndef NDEBUG
2413                              if (fabs(movement) >= 1.0e30)
2414                                   resetFakeBounds(-1000 - iSequence);
2415#endif
2416#ifdef CLP_DEBUG
2417                              if ((handler_->logLevel() & 32))
2418                                   printf("%d %d, new dj %g, alpha %g, movement %g\n",
2419                                          1, iSequence, value, alphaI, movement);
2420#endif
2421                              changeObj += movement * cost[iSequence];
2422                              matrix_->add(this, outputArray, iSequence, movement);
2423                         }
2424                    } else if (status == atUpperBound) {
2425                         double value = reducedCost[iSequence] - theta * alphaI;
2426                         reducedCost[iSequence] = value;
2427                         double movement = 0.0;
2428
2429                         if (value > tolerance) {
2430                              // to lower bound (if swap)
2431                              which[numberInfeasibilities++] = iSequence;
2432                              movement = lower[iSequence] - upper[iSequence];
2433#ifndef NDEBUG
2434                              if (fabs(movement) >= 1.0e30)
2435                                   resetFakeBounds(-1000 - iSequence);
2436#endif
2437#ifdef CLP_DEBUG
2438                              if ((handler_->logLevel() & 32))
2439                                   printf("%d %d, new dj %g, alpha %g, movement %g\n",
2440                                          1, iSequence, value, alphaI, movement);
2441#endif
2442                              changeObj += movement * cost[iSequence];
2443                              matrix_->add(this, outputArray, iSequence, movement);
2444                         }
2445                    } else if (status == isFree) {
2446                         double value = reducedCost[iSequence] - theta * alphaI;
2447                         reducedCost[iSequence] = value;
2448                    }
2449               }
2450          }
2451     } else  {
2452          double * COIN_RESTRICT solution = solutionRegion(0);
2453          double * COIN_RESTRICT reducedCost = djRegion(0);
2454          double * COIN_RESTRICT lower = lowerRegion(0);
2455          double * COIN_RESTRICT upper = upperRegion(0);
2456          const double * COIN_RESTRICT cost = costRegion(0);
2457          int * COIN_RESTRICT which;
2458          which = rowArray->getIndices();
2459          int iSequence;
2460          for (iSequence = 0; iSequence < numberRows_; iSequence++) {
2461               double value = reducedCost[iSequence];
2462
2463               Status status = getStatus(iSequence + numberColumns_);
2464               // more likely to be at upper bound ?
2465               if (status == atUpperBound) {
2466                    double movement = 0.0;
2467                    //#define NO_SWAP7
2468                    if (value > tolerance) {
2469                         // to lower bound (if swap)
2470                         // put back alpha
2471                         which[numberInfeasibilities++] = iSequence;
2472                         movement = lower[iSequence] - upper[iSequence];
2473#define TRY_SET_FAKE
2474#ifdef TRY_SET_FAKE
2475                         if (fabs(movement) > dualBound_) {
2476                           FakeBound bound = getFakeBound(iSequence + numberColumns_);
2477                           if (bound == ClpSimplexDual::noFake) {
2478                             setFakeBound(iSequence + numberColumns_, 
2479                                          ClpSimplexDual::lowerFake);
2480                             lower[iSequence] = upper[iSequence] - dualBound_;
2481                             assert (fabs(lower[iSequence])<1.0e30);
2482                             movement = lower[iSequence] - upper[iSequence];
2483                             numberFake_++;
2484#ifndef NDEBUG
2485                           } else {
2486                             if (fabs(movement) >= 1.0e30)
2487                               resetFakeBounds(-1000 - iSequence);
2488#endif
2489                           }
2490                         }
2491#endif
2492                         changeObj += movement * cost[iSequence];
2493                         outputArray->quickAdd(iSequence, -movement);
2494#ifndef NO_SWAP7
2495                    } else if (value > -tolerance) {
2496                         // at correct bound but may swap
2497                         FakeBound bound = getFakeBound(iSequence + numberColumns_);
2498                         if (bound == ClpSimplexDual::upperFake) {
2499                              movement = lower[iSequence] - upper[iSequence];
2500#ifndef NDEBUG
2501                              if (fabs(movement) >= 1.0e30)
2502                                resetFakeBounds(-1000 - iSequence);
2503#endif
2504                              setStatus(iSequence + numberColumns_, atLowerBound);
2505                              solution[iSequence] = lower[iSequence];
2506                              changeObj += movement * cost[iSequence];
2507                              //numberFake_--;
2508                              //setFakeBound(iSequence+numberColumns_,noFake);
2509                         }
2510#endif
2511                    }
2512               } else if (status == atLowerBound) {
2513                    double movement = 0.0;
2514
2515                    if (value < -tolerance) {
2516                         // to upper bound
2517                         // put back alpha
2518                         which[numberInfeasibilities++] = iSequence;
2519                         movement = upper[iSequence] - lower[iSequence];
2520#ifdef TRY_SET_FAKE
2521                         if (fabs(movement) > dualBound_) {
2522                           FakeBound bound = getFakeBound(iSequence + numberColumns_);
2523                           if (bound == ClpSimplexDual::noFake) {
2524                             setFakeBound(iSequence + numberColumns_, 
2525                                          ClpSimplexDual::upperFake);
2526                             upper[iSequence] = lower[iSequence] + dualBound_;
2527                             assert (fabs(upper[iSequence])<1.0e30);
2528                             movement = upper[iSequence] - lower[iSequence];
2529                             numberFake_++;
2530#ifndef NDEBUG
2531                           } else {
2532                             if (fabs(movement) >= 1.0e30)
2533                               resetFakeBounds(-1000 - iSequence);
2534#endif
2535                           }
2536                         }
2537#endif
2538                         changeObj += movement * cost[iSequence];
2539                         outputArray->quickAdd(iSequence, -movement);
2540#ifndef NO_SWAP7
2541                    } else if (value < tolerance) {
2542                         // at correct bound but may swap
2543                         FakeBound bound = getFakeBound(iSequence + numberColumns_);
2544                         if (bound == ClpSimplexDual::lowerFake) {
2545                              movement = upper[iSequence] - lower[iSequence];
2546#ifndef NDEBUG
2547                              if (fabs(movement) >= 1.0e30)
2548                                resetFakeBounds(-1000 - iSequence);
2549#endif
2550                              setStatus(iSequence + numberColumns_, atUpperBound);
2551                              solution[iSequence] = upper[iSequence];
2552                              changeObj += movement * cost[iSequence];
2553                              //numberFake_--;
2554                              //setFakeBound(iSequence+numberColumns_,noFake);
2555                         }
2556#endif
2557                    }
2558               }
2559          }
2560          // Do columns
2561          solution = solutionRegion(1);
2562          reducedCost = djRegion(1);
2563          lower = lowerRegion(1);
2564          upper = upperRegion(1);
2565          cost = costRegion(1);
2566          // set number of infeasibilities in row array
2567          numberRowInfeasibilities = numberInfeasibilities;
2568          rowArray->setNumElements(numberInfeasibilities);
2569          numberInfeasibilities = 0;
2570          which = columnArray->getIndices();
2571          for (iSequence = 0; iSequence < numberColumns_; iSequence++) {
2572               double value = reducedCost[iSequence];
2573
2574               Status status = getStatus(iSequence);
2575               if (status == atLowerBound) {
2576                    double movement = 0.0;
2577
2578                    if (value < -tolerance) {
2579                         // to upper bound
2580                         // put back alpha
2581                         which[numberInfeasibilities++] = iSequence;
2582                         movement = upper[iSequence] - lower[iSequence];
2583#ifdef TRY_SET_FAKE
2584                         if (fabs(movement) > dualBound_) {
2585                           FakeBound bound = getFakeBound(iSequence);
2586                           if (bound == ClpSimplexDual::noFake) {
2587                             setFakeBound(iSequence, 
2588                                          ClpSimplexDual::upperFake);
2589                             upper[iSequence] = lower[iSequence] + dualBound_;
2590                             assert (fabs(upper[iSequence])<1.0e30);
2591                             movement = upper[iSequence] - lower[iSequence];
2592                             numberFake_++;
2593#ifndef NDEBUG
2594                           } else {
2595                             if (fabs(movement) >= 1.0e30)
2596                               resetFakeBounds(-1000 - iSequence);
2597#endif
2598                           }
2599                         }
2600#endif
2601                         changeObj += movement * cost[iSequence];
2602                         matrix_->add(this, outputArray, iSequence, movement);
2603#ifndef NO_SWAP7
2604                    } else if (value < tolerance) {
2605                         // at correct bound but may swap
2606                         FakeBound bound = getFakeBound(iSequence);
2607                         if (bound == ClpSimplexDual::lowerFake) {
2608                              movement = upper[iSequence] - lower[iSequence];
2609#ifndef NDEBUG
2610                              if (fabs(movement) >= 1.0e30)
2611                                resetFakeBounds(-1000 - iSequence);
2612#endif
2613                              setStatus(iSequence, atUpperBound);
2614                              solution[iSequence] = upper[iSequence];
2615                              changeObj += movement * cost[iSequence];
2616                              //numberFake_--;
2617                              //setFakeBound(iSequence,noFake);
2618                         }
2619#endif
2620                    }
2621               } else if (status == atUpperBound) {
2622                    double movement = 0.0;
2623
2624                    if (value > tolerance) {
2625                         // to lower bound (if swap)
2626                         // put back alpha
2627                         which[numberInfeasibilities++] = iSequence;
2628                         movement = lower[iSequence] - upper[iSequence];
2629#ifdef TRY_SET_FAKE
2630                         if (fabs(movement) > dualBound_) {
2631                           FakeBound bound = getFakeBound(iSequence);
2632                           if (bound == ClpSimplexDual::noFake) {
2633                             setFakeBound(iSequence, 
2634                                          ClpSimplexDual::lowerFake);
2635                             lower[iSequence] = upper[iSequence] - dualBound_;
2636                             assert (fabs(lower[iSequence])<1.0e30);
2637                             movement = lower[iSequence] - upper[iSequence];
2638                             numberFake_++;
2639#ifndef NDEBUG
2640                           } else {
2641                             if (fabs(movement) >= 1.0e30)
2642                               resetFakeBounds(-1000 - iSequence);
2643#endif
2644                           }
2645                         }
2646#endif
2647                         changeObj += movement * cost[iSequence];
2648                         matrix_->add(this, outputArray, iSequence, movement);
2649#ifndef NO_SWAP7
2650                    } else if (value > -tolerance) {
2651                         // at correct bound but may swap
2652                         FakeBound bound = getFakeBound(iSequence);
2653                         if (bound == ClpSimplexDual::upperFake) {
2654                              movement = lower[iSequence] - upper[iSequence];
2655#ifndef NDEBUG
2656                              if (fabs(movement) >= 1.0e30)
2657                                resetFakeBounds(-1000 - iSequence);
2658#endif
2659                              setStatus(iSequence, atLowerBound);
2660                              solution[iSequence] = lower[iSequence];
2661                              changeObj += movement * cost[iSequence];
2662                              //numberFake_--;
2663                              //setFakeBound(iSequence,noFake);
2664                         }
2665#endif
2666                    }
2667               }
2668          }
2669     }
2670
2671#ifdef CLP_DEBUG
2672     if (fullRecompute && numberFake_ && (handler_->logLevel() & 16) != 0)
2673          printf("%d fake after full update\n", numberFake_);
2674#endif
2675     // set number of infeasibilities
2676     columnArray->setNumElements(numberInfeasibilities);
2677     numberInfeasibilities += numberRowInfeasibilities;
2678     if (fullRecompute) {
2679          // do actual flips
2680          flipBounds(rowArray, columnArray);
2681     }
2682     objectiveChange += changeObj;
2683     return numberInfeasibilities;
2684}
2685void
2686ClpSimplexDual::updateDualsInValuesPass(CoinIndexedVector * rowArray,
2687                                        CoinIndexedVector * columnArray,
2688                                        double theta)
2689{
2690
2691     // use a tighter tolerance except for all being okay
2692     double tolerance = dualTolerance_;
2693
2694     // Coding is very similar but we can save a bit by splitting
2695     // Do rows
2696     {
2697          int i;
2698          double * reducedCost = djRegion(0);
2699          double * work;
2700          int number;
2701          int * which;
2702          work = rowArray->denseVector();
2703          number = rowArray->getNumElements();
2704          which = rowArray->getIndices();
2705          for (i = 0; i < number; i++) {
2706               int iSequence = which[i];
2707               double alphaI = work[i];
2708               double value = reducedCost[iSequence] - theta * alphaI;
2709               work[i] = 0.0;
2710               reducedCost[iSequence] = value;
2711
2712               Status status = getStatus(iSequence + numberColumns_);
2713               // more likely to be at upper bound ?
2714               if (status == atUpperBound) {
2715
2716                    if (value > tolerance)
2717                         reducedCost[iSequence] = 0.0;
2718               } else if (status == atLowerBound) {
2719
2720                    if (value < -tolerance) {
2721                         reducedCost[iSequence] = 0.0;
2722                    }
2723               }
2724          }
2725     }
2726     rowArray->setNumElements(0);
2727
2728     // Do columns
2729     {
2730          int i;
2731          double * reducedCost = djRegion(1);
2732          double * work;
2733          int number;
2734          int * which;
2735          work = columnArray->denseVector();
2736          number = columnArray->getNumElements();
2737          which = columnArray->getIndices();
2738
2739          for (i = 0; i < number; i++) {
2740               int iSequence = which[i];
2741               double alphaI = work[i];
2742               double value = reducedCost[iSequence] - theta * alphaI;
2743               work[i] = 0.0;
2744               reducedCost[iSequence] = value;
2745
2746               Status status = getStatus(iSequence);
2747               if (status == atLowerBound) {
2748                    if (value < -tolerance)
2749                         reducedCost[iSequence] = 0.0;
2750               } else if (status == atUpperBound) {
2751                    if (value > tolerance)
2752                         reducedCost[iSequence] = 0.0;
2753               }
2754          }
2755     }
2756     columnArray->setNumElements(0);
2757}
2758/*
2759   Chooses dual pivot row
2760   Would be faster with separate region to scan
2761   and will have this (with square of infeasibility) when steepest
2762   For easy problems we can just choose one of the first rows we look at
2763*/
2764void
2765ClpSimplexDual::dualRow(int alreadyChosen)
2766{
2767     // get pivot row using whichever method it is
2768     int chosenRow = -1;
2769#ifdef FORCE_FOLLOW
2770     bool forceThis = false;
2771     if (!fpFollow && strlen(forceFile)) {
2772          fpFollow = fopen(forceFile, "r");
2773          assert (fpFollow);
2774     }
2775     if (fpFollow) {
2776          if (numberIterations_ <= force_iteration) {
2777               // read to next Clp0102
2778               char temp[300];
2779               while (fgets(temp, 250, fpFollow)) {
2780                    if (strncmp(temp, "Clp0102", 7))
2781                         continue;
2782                    char cin, cout;
2783                    sscanf(temp + 9, "%d%*f%*s%*c%c%d%*s%*c%c%d",
2784                           &force_iteration, &cin, &force_in, &cout, &force_out);
2785                    if (cin == 'R')
2786                         force_in += numberColumns_;
2787                    if (cout == 'R')
2788                         force_out += numberColumns_;
2789                    forceThis = true;
2790                    assert (numberIterations_ == force_iteration - 1);
2791                    printf("Iteration %d will force %d out and %d in\n",
2792                           force_iteration, force_out, force_in);
2793                    alreadyChosen = force_out;
2794                    break;
2795               }
2796          } else {
2797               // use old
2798               forceThis = true;
2799          }
2800          if (!forceThis) {
2801               fclose(fpFollow);
2802               fpFollow = NULL;
2803               forceFile = "";
2804          }
2805     }
2806#endif
2807     //double freeAlpha = 0.0;
2808     if (alreadyChosen < 0) {
2809          // first see if any free variables and put them in basis
2810          int nextFree = nextSuperBasic();
2811          //nextFree=-1; //off
2812          if (nextFree >= 0) {
2813               // unpack vector and find a good pivot
2814               unpack(rowArray_[1], nextFree);
2815               factorization_->updateColumn(rowArray_[2], rowArray_[1]);
2816
2817               double * work = rowArray_[1]->denseVector();
2818               int number = rowArray_[1]->getNumElements();
2819               int * which = rowArray_[1]->getIndices();
2820               double bestFeasibleAlpha = 0.0;
2821               int bestFeasibleRow = -1;
2822               double bestInfeasibleAlpha = 0.0;
2823               int bestInfeasibleRow = -1;
2824               int i;
2825
2826               for (i = 0; i < number; i++) {
2827                    int iRow = which[i];
2828                    double alpha = fabs(work[iRow]);
2829                    if (alpha > 1.0e-3) {
2830                         int iSequence = pivotVariable_[iRow];
2831                         double value = solution_[iSequence];
2832                         double lower = lower_[iSequence];
2833                         double upper = upper_[iSequence];
2834                         double infeasibility = 0.0;
2835                         if (value > upper)
2836                              infeasibility = value - upper;
2837                         else if (value < lower)
2838                              infeasibility = lower - value;
2839                         if (infeasibility * alpha > bestInfeasibleAlpha && alpha > 1.0e-1) {
2840                              if (!flagged(iSequence)) {
2841                                   bestInfeasibleAlpha = infeasibility * alpha;
2842                                   bestInfeasibleRow = iRow;
2843                              }
2844                         }
2845                         if (alpha > bestFeasibleAlpha && (lower > -1.0e20 || upper < 1.0e20)) {
2846                              bestFeasibleAlpha = alpha;
2847                              bestFeasibleRow = iRow;
2848                         }
2849                    }
2850               }
2851               if (bestInfeasibleRow >= 0)
2852                    chosenRow = bestInfeasibleRow;
2853               else if (bestFeasibleAlpha > 1.0e-2)
2854                    chosenRow = bestFeasibleRow;
2855               if (chosenRow >= 0) {
2856                    pivotRow_ = chosenRow;
2857                    //freeAlpha = work[chosenRow];
2858               }
2859               rowArray_[1]->clear();
2860          }
2861     } else {
2862          // in values pass
2863          chosenRow = alreadyChosen;
2864#ifdef FORCE_FOLLOW
2865          if(forceThis) {
2866               alreadyChosen = -1;
2867               chosenRow = -1;
2868               for (int i = 0; i < numberRows_; i++) {
2869                    if (pivotVariable_[i] == force_out) {
2870                         chosenRow = i;
2871                         break;
2872                    }
2873               }
2874               assert (chosenRow >= 0);
2875          }
2876#endif
2877          pivotRow_ = chosenRow;
2878     }
2879     if (chosenRow < 0)
2880          pivotRow_ = dualRowPivot_->pivotRow();
2881
2882     if (pivotRow_ >= 0) {
2883          sequenceOut_ = pivotVariable_[pivotRow_];
2884          valueOut_ = solution_[sequenceOut_];
2885          lowerOut_ = lower_[sequenceOut_];
2886          upperOut_ = upper_[sequenceOut_];
2887          if (alreadyChosen < 0) {
2888               // if we have problems we could try other way and hope we get a
2889               // zero pivot?
2890               if (valueOut_ > upperOut_) {
2891                    directionOut_ = -1;
2892                    dualOut_ = valueOut_ - upperOut_;
2893               } else if (valueOut_ < lowerOut_) {
2894                    directionOut_ = 1;
2895                    dualOut_ = lowerOut_ - valueOut_;
2896               } else {
2897#if 1
2898                    // odd (could be free) - it's feasible - go to nearest
2899                    if (valueOut_ - lowerOut_ < upperOut_ - valueOut_) {
2900                         directionOut_ = 1;
2901                         dualOut_ = lowerOut_ - valueOut_;
2902                    } else {
2903                         directionOut_ = -1;
2904                         dualOut_ = valueOut_ - upperOut_;
2905                    }
2906#else
2907                    // odd (could be free) - it's feasible - improve obj
2908                    printf("direction from alpha of %g is %d\n",
2909                           freeAlpha, freeAlpha > 0.0 ? 1 : -1);
2910                    if (valueOut_ - lowerOut_ > 1.0e20)
2911                         freeAlpha = 1.0;
2912                    else if(upperOut_ - valueOut_ > 1.0e20)
2913                         freeAlpha = -1.0;
2914                    //if (valueOut_-lowerOut_<upperOut_-valueOut_) {
2915                    if (freeAlpha < 0.0) {
2916                         directionOut_ = 1;
2917                         dualOut_ = lowerOut_ - valueOut_;
2918                    } else {
2919                         directionOut_ = -1;
2920                         dualOut_ = valueOut_ - upperOut_;
2921                    }
2922                    printf("direction taken %d - bounds %g %g %g\n",
2923                           directionOut_, lowerOut_, valueOut_, upperOut_);
2924#endif
2925               }
2926#ifdef CLP_DEBUG
2927               assert(dualOut_ >= 0.0);
2928#endif
2929          } else {
2930               // in values pass so just use sign of dj
2931               // We don't want to go through any barriers so set dualOut low
2932               // free variables will never be here
2933               dualOut_ = 1.0e-6;
2934               if (dj_[sequenceOut_] > 0.0) {
2935                    // this will give a -1 in pivot row (as slacks are -1.0)
2936                    directionOut_ = 1;
2937               } else {
2938                    directionOut_ = -1;
2939               }
2940          }
2941     }
2942     return ;
2943}
2944// Checks if any fake bounds active - if so returns number and modifies
2945// dualBound_ and everything.
2946// Free variables will be left as free
2947// Returns number of bounds changed if >=0
2948// Returns -1 if not initialize and no effect
2949// Fills in changeVector which can be used to see if unbounded
2950// and cost of change vector
2951int
2952ClpSimplexDual::changeBounds(int initialize,
2953                             CoinIndexedVector * outputArray,
2954                             double & changeCost)
2955{
2956     numberFake_ = 0;
2957     if (!initialize) {
2958          int numberInfeasibilities;
2959          double newBound;
2960          newBound = 5.0 * dualBound_;
2961          numberInfeasibilities = 0;
2962          changeCost = 0.0;
2963          // put back original bounds and then check
2964          createRim1(false);
2965          int iSequence;
2966          // bounds will get bigger - just look at ones at bounds
2967          for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
2968               double lowerValue = lower_[iSequence];
2969               double upperValue = upper_[iSequence];
2970               double value = solution_[iSequence];
2971               setFakeBound(iSequence, ClpSimplexDual::noFake);
2972               switch(getStatus(iSequence)) {
2973
2974               case basic:
2975               case ClpSimplex::isFixed:
2976                    break;
2977               case isFree:
2978               case superBasic:
2979                    break;
2980               case atUpperBound:
2981                    if (fabs(value - upperValue) > primalTolerance_)
2982                         numberInfeasibilities++;
2983                    break;
2984               case atLowerBound:
2985                    if (fabs(value - lowerValue) > primalTolerance_)
2986                         numberInfeasibilities++;
2987                    break;
2988               }
2989          }
2990          // If dual infeasible then carry on
2991          if (numberInfeasibilities) {
2992               handler_->message(CLP_DUAL_CHECKB, messages_)
2993                         << newBound
2994                         << CoinMessageEol;
2995               int iSequence;
2996               for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
2997                    double lowerValue = lower_[iSequence];
2998                    double upperValue = upper_[iSequence];
2999                    double newLowerValue;
3000                    double newUpperValue;
3001                    Status status = getStatus(iSequence);
3002                    if (status == atUpperBound ||
3003                              status == atLowerBound) {
3004                         double value = solution_[iSequence];
3005                         if (value - lowerValue <= upperValue - value) {
3006                              newLowerValue = CoinMax(lowerValue, value - 0.666667 * newBound);
3007                              newUpperValue = CoinMin(upperValue, newLowerValue + newBound);
3008                         } else {
3009                              newUpperValue = CoinMin(upperValue, value + 0.666667 * newBound);
3010                              newLowerValue = CoinMax(lowerValue, newUpperValue - newBound);
3011                         }
3012                         lower_[iSequence] = newLowerValue;
3013                         upper_[iSequence] = newUpperValue;
3014                         if (newLowerValue > lowerValue) {
3015                              if (newUpperValue < upperValue) {
3016                                   setFakeBound(iSequence, ClpSimplexDual::bothFake);
3017#ifdef CLP_INVESTIGATE
3018                                   abort(); // No idea what should happen here - I have never got here
3019#endif
3020                                   numberFake_++;
3021                              } else {
3022                                   setFakeBound(iSequence, ClpSimplexDual::lowerFake);
3023                                   numberFake_++;
3024                              }
3025                         } else {
3026                              if (newUpperValue < upperValue) {
3027                                   setFakeBound(iSequence, ClpSimplexDual::upperFake);
3028                                   numberFake_++;
3029                              }
3030                         }
3031                         if (status == atUpperBound)
3032                              solution_[iSequence] = newUpperValue;
3033                         else
3034                              solution_[iSequence] = newLowerValue;
3035                         double movement = solution_[iSequence] - value;
3036                         if (movement && outputArray) {
3037                              if (iSequence >= numberColumns_) {
3038                                   outputArray->quickAdd(iSequence, -movement);
3039                                   changeCost += movement * cost_[iSequence];
3040                              } else {
3041                                   matrix_->add(this, outputArray, iSequence, movement);
3042                                   changeCost += movement * cost_[iSequence];
3043                              }
3044                         }
3045                    }
3046               }
3047               dualBound_ = newBound;
3048          } else {
3049               numberInfeasibilities = -1;
3050          }
3051          return numberInfeasibilities;
3052     } else if (initialize == 1 || initialize == 3) {
3053          int iSequence;
3054          if (initialize == 3) {
3055               for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
3056                    setFakeBound(iSequence, ClpSimplexDual::noFake);
3057               }
3058          }
3059          double testBound = 0.999999 * dualBound_;
3060          for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
3061               Status status = getStatus(iSequence);
3062               if (status == atUpperBound ||
3063                         status == atLowerBound) {
3064                    double lowerValue = lower_[iSequence];
3065                    double upperValue = upper_[iSequence];
3066                    double value = solution_[iSequence];
3067                    if (lowerValue > -largeValue_ || upperValue < largeValue_) {
3068                         if (true || lowerValue - value > -0.5 * dualBound_ ||
3069                                   upperValue - value < 0.5 * dualBound_) {
3070                              if (fabs(lowerValue - value) <= fabs(upperValue - value)) {
3071                                   if (upperValue > lowerValue + testBound) {
3072                                        if (getFakeBound(iSequence) == ClpSimplexDual::noFake)
3073                                             numberFake_++;
3074                                        upper_[iSequence] = lowerValue + dualBound_;
3075                                        setFakeBound(iSequence, ClpSimplexDual::upperFake);
3076                                   }
3077                              } else {
3078                                   if (lowerValue < upperValue - testBound) {
3079                                        if (getFakeBound(iSequence) == ClpSimplexDual::noFake)
3080                                             numberFake_++;
3081                                        lower_[iSequence] = upperValue - dualBound_;
3082                                        setFakeBound(iSequence, ClpSimplexDual::lowerFake);
3083                                   }
3084                              }
3085                         } else {
3086                              if (getFakeBound(iSequence) == ClpSimplexDual::noFake)
3087                                   numberFake_++;
3088                              lower_[iSequence] = -0.5 * dualBound_;
3089                              upper_[iSequence] = 0.5 * dualBound_;
3090                              setFakeBound(iSequence, ClpSimplexDual::bothFake);
3091                              abort();
3092                         }
3093                         if (status == atUpperBound)
3094                              solution_[iSequence] = upper_[iSequence];
3095                         else
3096                              solution_[iSequence] = lower_[iSequence];
3097                    } else {
3098                         // set non basic free variables to fake bounds
3099                         // I don't think we should ever get here
3100                         CoinAssert(!("should not be here"));
3101                         lower_[iSequence] = -0.5 * dualBound_;
3102                         upper_[iSequence] = 0.5 * dualBound_;
3103                         setFakeBound(iSequence, ClpSimplexDual::bothFake);
3104                         numberFake_++;
3105                         setStatus(iSequence, atUpperBound);
3106                         solution_[iSequence] = 0.5 * dualBound_;
3107                    }
3108               } else if (status == basic) {
3109                    // make sure not at fake bound and bounds correct
3110                    setFakeBound(iSequence, ClpSimplexDual::noFake);
3111                    double gap = upper_[iSequence] - lower_[iSequence];
3112                    if (gap > 0.5 * dualBound_ && gap < 2.0 * dualBound_) {
3113                         if (iSequence < numberColumns_) {
3114                              if (columnScale_) {
3115                                   double multiplier = rhsScale_ * inverseColumnScale_[iSequence];
3116                                   // lower
3117                                   double value = columnLower_[iSequence];
3118                                   if (value > -1.0e30) {
3119                                        value *= multiplier;
3120                                   }
3121                                   lower_[iSequence] = value;
3122                                   // upper
3123                                   value = columnUpper_[iSequence];
3124                                   if (value < 1.0e30) {
3125                                        value *= multiplier;
3126                                   }
3127                                   upper_[iSequence] = value;
3128                              } else {
3129                                   lower_[iSequence] = columnLower_[iSequence];;
3130                                   upper_[iSequence] = columnUpper_[iSequence];;
3131                              }
3132                         } else {
3133                              int iRow = iSequence - numberColumns_;
3134                              if (rowScale_) {
3135                                   // lower
3136                                   double multiplier = rhsScale_ * rowScale_[iRow];
3137                                   double value = rowLower_[iRow];
3138                                   if (value > -1.0e30) {
3139                                        value *= multiplier;
3140                                   }
3141                                   lower_[iSequence] = value;
3142                                   // upper
3143                                   value = rowUpper_[iRow];
3144                                   if (value < 1.0e30) {
3145                                        value *= multiplier;
3146                                   }
3147                                   upper_[iSequence] = value;
3148                              } else {
3149                                   lower_[iSequence] = rowLower_[iRow];;
3150                                   upper_[iSequence] = rowUpper_[iRow];;
3151                              }
3152                         }
3153                    }
3154               }
3155          }
3156
3157          return 1;
3158     } else {
3159          // just reset changed ones
3160          if (columnScale_) {
3161               int iSequence;
3162               for (iSequence = 0; iSequence < numberColumns_; iSequence++) {
3163                    FakeBound fakeStatus = getFakeBound(iSequence);
3164                    if (fakeStatus != noFake) {
3165                         if ((static_cast<int> (fakeStatus) & 1) != 0) {
3166                              // lower
3167                              double value = columnLower_[iSequence];
3168                              if (value > -1.0e30) {
3169                                   double multiplier = rhsScale_ * inverseColumnScale_[iSequence];
3170                                   value *= multiplier;
3171                              }
3172                              columnLowerWork_[iSequence] = value;
3173                         }
3174                         if ((static_cast<int> (fakeStatus) & 2) != 0) {
3175                              // upper
3176                              double value = columnUpper_[iSequence];
3177                              if (value < 1.0e30) {
3178                                   double multiplier = rhsScale_ * inverseColumnScale_[iSequence];
3179                                   value *= multiplier;
3180                              }
3181                              columnUpperWork_[iSequence] = value;
3182                         }
3183                    }
3184               }
3185               for (iSequence = 0; iSequence < numberRows_; iSequence++) {
3186                    FakeBound fakeStatus = getFakeBound(iSequence + numberColumns_);
3187                    if (fakeStatus != noFake) {
3188                         if ((static_cast<int> (fakeStatus) & 1) != 0) {
3189                              // lower
3190                              double value = rowLower_[iSequence];
3191                              if (value > -1.0e30) {
3192                                   double multiplier = rhsScale_ * rowScale_[iSequence];
3193                                   value *= multiplier;
3194                              }
3195                              rowLowerWork_[iSequence] = value;
3196                         }
3197                         if ((static_cast<int> (fakeStatus) & 2) != 0) {
3198                              // upper
3199                              double value = rowUpper_[iSequence];
3200                              if (value < 1.0e30) {
3201                                   double multiplier = rhsScale_ * rowScale_[iSequence];
3202                                   value *= multiplier;
3203                              }
3204                              rowUpperWork_[iSequence] = value;
3205                         }
3206                    }
3207               }
3208          } else {
3209               int iSequence;
3210               for (iSequence = 0; iSequence < numberColumns_; iSequence++) {
3211                    FakeBound fakeStatus = getFakeBound(iSequence);
3212                    if ((static_cast<int> (fakeStatus) & 1) != 0) {
3213                         // lower
3214                         columnLowerWork_[iSequence] = columnLower_[iSequence];
3215                    }
3216                    if ((static_cast<int> (fakeStatus) & 2) != 0) {
3217                         // upper
3218                         columnUpperWork_[iSequence] = columnUpper_[iSequence];
3219                    }
3220               }
3221               for (iSequence = 0; iSequence < numberRows_; iSequence++) {
3222                    FakeBound fakeStatus = getFakeBound(iSequence + numberColumns_);
3223                    if ((static_cast<int> (fakeStatus) & 1) != 0) {
3224                         // lower
3225                         rowLowerWork_[iSequence] = rowLower_[iSequence];
3226                    }
3227                    if ((static_cast<int> (fakeStatus) & 2) != 0) {
3228                         // upper
3229                         rowUpperWork_[iSequence] = rowUpper_[iSequence];
3230                    }
3231               }
3232          }
3233          return 0;
3234     }
3235}
3236int
3237ClpSimplexDual::dualColumn0(const CoinIndexedVector * rowArray,
3238                            const CoinIndexedVector * columnArray,
3239                            CoinIndexedVector * spareArray,
3240                            double acceptablePivot,
3241                            double & upperReturn, double &bestReturn, double & badFree)
3242{
3243     // do first pass to get possibles
3244     double * spare = spareArray->denseVector();
3245     int * index = spareArray->getIndices();
3246     const double * work;
3247     int number;
3248     const int * which;
3249     const double * reducedCost;
3250     // We can also see if infeasible or pivoting on free
3251     double tentativeTheta = 1.0e15;
3252     double upperTheta = 1.0e31;
3253     double freePivot = acceptablePivot;
3254     double bestPossible = 0.0;
3255     int numberRemaining = 0;
3256     int i;
3257     badFree = 0.0;
3258     if ((moreSpecialOptions_ & 8) != 0) {
3259          // No free or super basic
3260          double multiplier[] = { -1.0, 1.0};
3261          double dualT = - dualTolerance_;
3262          for (int iSection = 0; iSection < 2; iSection++) {
3263
3264               int addSequence;
3265               unsigned char * statusArray;
3266               if (!iSection) {
3267                    work = rowArray->denseVector();
3268                    number = rowArray->getNumElements();
3269                    which = rowArray->getIndices();
3270                    reducedCost = rowReducedCost_;
3271                    addSequence = numberColumns_;
3272                    statusArray = status_ + numberColumns_;
3273               } else {
3274                    work = columnArray->denseVector();
3275                    number = columnArray->getNumElements();
3276                    which = columnArray->getIndices();
3277                    reducedCost = reducedCostWork_;
3278                    addSequence = 0;
3279                    statusArray = status_;
3280               }
3281
3282               for (i = 0; i < number; i++) {
3283                    int iSequence = which[i];
3284                    double alpha;
3285                    double oldValue;
3286                    double value;
3287
3288                    assert (getStatus(iSequence + addSequence) != isFree
3289                            && getStatus(iSequence + addSequence) != superBasic);
3290                    int iStatus = (statusArray[iSequence] & 3) - 1;
3291                    if (iStatus) {
3292                         double mult = multiplier[iStatus-1];
3293                         alpha = work[i] * mult;
3294                         if (alpha > 0.0) {
3295                              oldValue = reducedCost[iSequence] * mult;
3296                              value = oldValue - tentativeTheta * alpha;
3297                              if (value < dualT) {
3298                                   bestPossible = CoinMax(bestPossible, alpha);
3299                                   value = oldValue - upperTheta * alpha;
3300                                   if (value < dualT && alpha >= acceptablePivot) {
3301                                        upperTheta = (oldValue - dualT) / alpha;
3302                                        //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
3303                                   }
3304                                   // add to list
3305                                   spare[numberRemaining] = alpha * mult;
3306                                   index[numberRemaining++] = iSequence + addSequence;
3307                              }
3308                         }
3309                    }
3310               }
3311          }
3312     } else {
3313          // some free or super basic
3314          for (int iSection = 0; iSection < 2; iSection++) {
3315
3316               int addSequence;
3317
3318               if (!iSection) {
3319                    work = rowArray->denseVector();
3320                    number = rowArray->getNumElements();
3321                    which = rowArray->getIndices();
3322                    reducedCost = rowReducedCost_;
3323                    addSequence = numberColumns_;
3324               } else {
3325                    work = columnArray->denseVector();
3326                    number = columnArray->getNumElements();
3327                    which = columnArray->getIndices();
3328                    reducedCost = reducedCostWork_;
3329                    addSequence = 0;
3330               }
3331
3332               for (i = 0; i < number; i++) {
3333                    int iSequence = which[i];
3334                    double alpha;
3335                    double oldValue;
3336                    double value;
3337                    bool keep;
3338
3339                    switch(getStatus(iSequence + addSequence)) {
3340
3341                    case basic:
3342                    case ClpSimplex::isFixed:
3343                         break;
3344                    case isFree:
3345                    case superBasic:
3346                         alpha = work[i];
3347                         bestPossible = CoinMax(bestPossible, fabs(alpha));
3348                         oldValue = reducedCost[iSequence];
3349                         // If free has to be very large - should come in via dualRow
3350                         //if (getStatus(iSequence+addSequence)==isFree&&fabs(alpha)<1.0e-3)
3351                         //break;
3352                         if (oldValue > dualTolerance_) {
3353                              keep = true;
3354                         } else if (oldValue < -dualTolerance_) {
3355                              keep = true;
3356                         } else {
3357                              if (fabs(alpha) > CoinMax(10.0 * acceptablePivot, 1.0e-5)) {
3358                                   keep = true;
3359                              } else {
3360                                   keep = false;
3361                                   badFree = CoinMax(badFree, fabs(alpha));
3362                              }
3363                         }
3364                         if (keep) {
3365                              // free - choose largest
3366                              if (fabs(alpha) > freePivot) {
3367                                   freePivot = fabs(alpha);
3368                                   sequenceIn_ = iSequence + addSequence;
3369                                   theta_ = oldValue / alpha;
3370                                   alpha_ = alpha;
3371                              }
3372                         }
3373                         break;
3374                    case atUpperBound:
3375                         alpha = work[i];
3376                         oldValue = reducedCost[iSequence];
3377                         value = oldValue - tentativeTheta * alpha;
3378                         //assert (oldValue<=dualTolerance_*1.0001);
3379                         if (value > dualTolerance_) {
3380                              bestPossible = CoinMax(bestPossible, -alpha);
3381                              value = oldValue - upperTheta * alpha;
3382                              if (value > dualTolerance_ && -alpha >= acceptablePivot) {
3383                                   upperTheta = (oldValue - dualTolerance_) / alpha;
3384                                   //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
3385                              }
3386                              // add to list
3387                              spare[numberRemaining] = alpha;
3388                              index[numberRemaining++] = iSequence + addSequence;
3389                         }
3390                         break;
3391                    case atLowerBound:
3392                         alpha = work[i];
3393                         oldValue = reducedCost[iSequence];
3394                         value = oldValue - tentativeTheta * alpha;
3395                         //assert (oldValue>=-dualTolerance_*1.0001);
3396                         if (value < -dualTolerance_) {
3397                              bestPossible = CoinMax(bestPossible, alpha);
3398                              value = oldValue - upperTheta * alpha;
3399                              if (value < -dualTolerance_ && alpha >= acceptablePivot) {
3400                                   upperTheta = (oldValue + dualTolerance_) / alpha;
3401                                   //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
3402                              }
3403                              // add to list
3404                              spare[numberRemaining] = alpha;
3405                              index[numberRemaining++] = iSequence + addSequence;
3406                         }
3407                         break;
3408                    }
3409               }
3410          }
3411     }
3412     upperReturn = upperTheta;
3413     bestReturn = bestPossible;
3414     return numberRemaining;
3415}
3416/*
3417   Row array has row part of pivot row (as duals so sign may be switched)
3418   Column array has column part.
3419   This chooses pivot column.
3420   Spare array will be needed when we start getting clever.
3421   We will check for basic so spare array will never overflow.
3422   If necessary will modify costs
3423*/
3424double
3425ClpSimplexDual::dualColumn(CoinIndexedVector * rowArray,
3426                           CoinIndexedVector * columnArray,
3427                           CoinIndexedVector * spareArray,
3428                           CoinIndexedVector * spareArray2,
3429                           double acceptablePivot,
3430                           CoinBigIndex * /*dubiousWeights*/)
3431{
3432     int numberPossiblySwapped = 0;
3433     int numberRemaining = 0;
3434
3435     double totalThru = 0.0; // for when variables flip
3436     //double saveAcceptable=acceptablePivot;
3437     //acceptablePivot=1.0e-9;
3438
3439     double bestEverPivot = acceptablePivot;
3440     int lastSequence = -1;
3441     double lastPivot = 0.0;
3442     double upperTheta;
3443     double newTolerance = dualTolerance_;
3444     //newTolerance = dualTolerance_+1.0e-6*dblParam_[ClpDualTolerance];
3445     // will we need to increase tolerance
3446     //bool thisIncrease = false;
3447     // If we think we need to modify costs (not if something from broad sweep)
3448     bool modifyCosts = false;
3449     // Increase in objective due to swapping bounds (may be negative)
3450     double increaseInObjective = 0.0;
3451
3452     // use spareArrays to put ones looked at in
3453     // we are going to flip flop between
3454     int iFlip = 0;
3455     // Possible list of pivots
3456     int interesting[2];
3457     // where possible swapped ones are
3458     int swapped[2];
3459     // for zeroing out arrays after
3460     int marker[2][2];
3461     // pivot elements
3462     double * array[2], * spare, * spare2;
3463     // indices
3464     int * indices[2], * index, * index2;
3465     spareArray2->clear();
3466     array[0] = spareArray->denseVector();
3467     indices[0] = spareArray->getIndices();
3468     spare = array[0];
3469     index = indices[0];
3470     array[1] = spareArray2->denseVector();
3471     indices[1] = spareArray2->getIndices();
3472     int i;
3473
3474     // initialize lists
3475     for (i = 0; i < 2; i++) {
3476          interesting[i] = 0;
3477          swapped[i] = numberColumns_;
3478          marker[i][0] = 0;
3479          marker[i][1] = numberColumns_;
3480     }
3481     /*
3482       First we get a list of possible pivots.  We can also see if the
3483       problem looks infeasible or whether we want to pivot in free variable.
3484       This may make objective go backwards but can only happen a finite
3485       number of times and I do want free variables basic.
3486
3487       Then we flip back and forth.  At the start of each iteration
3488       interesting[iFlip] should have possible candidates and swapped[iFlip]
3489       will have pivots if we decide to take a previous pivot.
3490       At end of each iteration interesting[1-iFlip] should have
3491       candidates if we go through this theta and swapped[1-iFlip]
3492       pivots if we don't go through.
3493
3494       At first we increase theta and see what happens.  We start
3495       theta at a reasonable guess.  If in right area then we do bit by bit.
3496
3497      */
3498
3499     // do first pass to get possibles
3500     upperTheta = 1.0e31;
3501     double bestPossible = 0.0;
3502     double badFree = 0.0;
3503     alpha_ = 0.0;
3504     if (spareIntArray_[0] >= 0) {
3505          numberRemaining = dualColumn0(rowArray, columnArray, spareArray,
3506                                        acceptablePivot, upperTheta, bestPossible, badFree);
3507     } else {
3508          // already done
3509          numberRemaining = spareArray->getNumElements();
3510          spareArray->setNumElements(0);
3511          upperTheta = spareDoubleArray_[0];
3512          bestPossible = spareDoubleArray_[1];
3513          if (spareIntArray_[0] == -1) {
3514               theta_ = spareDoubleArray_[2];
3515               alpha_ = spareDoubleArray_[3];
3516               sequenceIn_ = spareIntArray_[1];
3517          } else {
3518#if 0
3519               int n = numberRemaining;
3520               double u = upperTheta;
3521               double b = bestPossible;
3522               upperTheta = 1.0e31;
3523               bestPossible = 0.0;
3524               numberRemaining = dualColumn0(rowArray, columnArray, spareArray,
3525                                             acceptablePivot, upperTheta, bestPossible, badFree);
3526               assert (n == numberRemaining);
3527               assert (fabs(b - bestPossible) < 1.0e-7);
3528               assert (fabs(u - upperTheta) < 1.0e-7);
3529#endif
3530          }
3531     }
3532     // switch off
3533     spareIntArray_[0] = 0;
3534     // We can also see if infeasible or pivoting on free
3535     double tentativeTheta = 1.0e25;
3536     interesting[0] = numberRemaining;
3537     marker[0][0] = numberRemaining;
3538
3539     if (!numberRemaining && sequenceIn_ < 0)
3540          return 0.0; // Looks infeasible
3541
3542     // If sum of bad small pivots too much
3543#define MORE_CAREFUL
3544#ifdef MORE_CAREFUL
3545     bool badSumPivots = false;
3546#endif
3547     if (sequenceIn_ >= 0) {
3548          // free variable - always choose
3549     } else {
3550
3551          theta_ = 1.0e50;
3552          // now flip flop between spare arrays until reasonable theta
3553          tentativeTheta = CoinMax(10.0 * upperTheta, 1.0e-7);
3554
3555          // loops increasing tentative theta until can't go through
3556
3557          while (tentativeTheta < 1.0e22) {
3558               double thruThis = 0.0;
3559
3560               double bestPivot = acceptablePivot;
3561               int bestSequence = -1;
3562
3563               numberPossiblySwapped = numberColumns_;
3564               numberRemaining = 0;
3565
3566               upperTheta = 1.0e50;
3567
3568               spare = array[iFlip];
3569               index = indices[iFlip];
3570               spare2 = array[1-iFlip];
3571               index2 = indices[1-iFlip];
3572
3573               // try 3 different ways
3574               // 1 bias increase by ones with slightly wrong djs
3575               // 2 bias by all
3576               // 3 bias by all - tolerance
3577#define TRYBIAS 3
3578
3579
3580               double increaseInThis = 0.0; //objective increase in this loop
3581
3582               for (i = 0; i < interesting[iFlip]; i++) {
3583                    int iSequence = index[i];
3584                    double alpha = spare[i];
3585                    double oldValue = dj_[iSequence];
3586                    double value = oldValue - tentativeTheta * alpha;
3587
3588                    if (alpha < 0.0) {
3589                         //at upper bound
3590                         if (value > newTolerance) {
3591                              double range = upper_[iSequence] - lower_[iSequence];
3592                              thruThis -= range * alpha;
3593#if TRYBIAS==1
3594                              if (oldValue > 0.0)
3595                                   increaseInThis -= oldValue * range;
3596#elif TRYBIAS==2
3597                              increaseInThis -= oldValue * range;
3598#else
3599                              increaseInThis -= (oldValue + dualTolerance_) * range;
3600#endif
3601                              // goes on swapped list (also means candidates if too many)
3602                              spare2[--numberPossiblySwapped] = alpha;
3603                              index2[numberPossiblySwapped] = iSequence;
3604                              if (fabs(alpha) > bestPivot) {
3605                                   bestPivot = fabs(alpha);
3606                                   bestSequence = numberPossiblySwapped;
3607                              }
3608                         } else {
3609                              value = oldValue - upperTheta * alpha;
3610                              if (value > newTolerance && -alpha >= acceptablePivot)
3611                                   upperTheta = (oldValue - newTolerance) / alpha;
3612                              spare2[numberRemaining] = alpha;
3613                              index2[numberRemaining++] = iSequence;
3614                         }
3615                    } else {
3616                         // at lower bound
3617                         if (value < -newTolerance) {
3618                              double range = upper_[iSequence] - lower_[iSequence];
3619                              thruThis += range * alpha;
3620                              //?? is this correct - and should we look at good ones
3621#if TRYBIAS==1
3622                              if (oldValue < 0.0)
3623                                   increaseInThis += oldValue * range;
3624#elif TRYBIAS==2
3625                              increaseInThis += oldValue * range;
3626#else
3627                              increaseInThis += (oldValue - dualTolerance_) * range;
3628#endif
3629                              // goes on swapped list (also means candidates if too many)
3630                              spare2[--numberPossiblySwapped] = alpha;
3631                              index2[numberPossiblySwapped] = iSequence;
3632                              if (fabs(alpha) > bestPivot) {
3633                                   bestPivot = fabs(alpha);
3634                                   bestSequence = numberPossiblySwapped;
3635                              }
3636                         } else {
3637                              value = oldValue - upperTheta * alpha;
3638                              if (value < -newTolerance && alpha >= acceptablePivot)
3639                                   upperTheta = (oldValue + newTolerance) / alpha;
3640                              spare2[numberRemaining] = alpha;
3641                              index2[numberRemaining++] = iSequence;
3642                         }
3643                    }
3644               }
3645               swapped[1-iFlip] = numberPossiblySwapped;
3646               interesting[1-iFlip] = numberRemaining;
3647               marker[1-iFlip][0] = CoinMax(marker[1-iFlip][0], numberRemaining);
3648               marker[1-iFlip][1] = CoinMin(marker[1-iFlip][1], numberPossiblySwapped);
3649
3650               double check = fabs(totalThru+thruThis);
3651               // add a bit
3652               check += 1.0e-8+1.0e-10*check;
3653               if (check >= fabs(dualOut_) ||
3654                         increaseInObjective + increaseInThis < 0.0) {
3655                    // We should be pivoting in this batch
3656                    // so compress down to this lot
3657                    numberRemaining = 0;
3658                    for (i = numberColumns_ - 1; i >= swapped[1-iFlip]; i--) {
3659                         spare[numberRemaining] = spare2[i];
3660                         index[numberRemaining++] = index2[i];
3661                    }
3662                    interesting[iFlip] = numberRemaining;
3663                    int iTry;
3664#define MAXTRY 100
3665                    // first get ratio with tolerance
3666                    for (iTry = 0; iTry < MAXTRY; iTry++) {
3667
3668                         upperTheta = 1.0e50;
3669                         numberPossiblySwapped = numberColumns_;
3670                         numberRemaining = 0;
3671
3672                         increaseInThis = 0.0; //objective increase in this loop
3673
3674                         thruThis = 0.0;
3675
3676                         spare = array[iFlip];
3677                         index = indices[iFlip];
3678                         spare2 = array[1-iFlip];
3679                         index2 = indices[1-iFlip];
3680                         for (i = 0; i < interesting[iFlip]; i++) {
3681                              int iSequence = index[i];
3682                              double alpha = spare[i];
3683                              double oldValue = dj_[iSequence];
3684                              double value = oldValue - upperTheta * alpha;
3685
3686                              if (alpha < 0.0) {
3687                                   //at upper bound
3688                                   if (value > newTolerance) {
3689                                        if (-alpha >= acceptablePivot) {
3690                                             upperTheta = (oldValue - newTolerance) / alpha;
3691                                        }
3692                                   }
3693                              } else {
3694                                   // at lower bound
3695                                   if (value < -newTolerance) {
3696                                        if (alpha >= acceptablePivot) {
3697                                             upperTheta = (oldValue + newTolerance) / alpha;
3698                                        }
3699                                   }
3700                              }
3701                         }
3702                         bestPivot = acceptablePivot;
3703                         sequenceIn_ = -1;
3704#ifdef DUBIOUS_WEIGHTS
3705                         double bestWeight = COIN_DBL_MAX;
3706#endif
3707                         double largestPivot = acceptablePivot;
3708                         // now choose largest and sum all ones which will go through
3709                         //printf("XX it %d number %d\n",numberIterations_,interesting[iFlip]);
3710                         // Sum of bad small pivots
3711#ifdef MORE_CAREFUL
3712                         double sumBadPivots = 0.0;
3713                         badSumPivots = false;
3714#endif
3715                         // Make sure upperTheta will work (-O2 and above gives problems)
3716                         upperTheta *= 1.0000000001;
3717                         for (i = 0; i < interesting[iFlip]; i++) {
3718                              int iSequence = index[i];
3719                              double alpha = spare[i];
3720                              double value = dj_[iSequence] - upperTheta * alpha;
3721                              double badDj = 0.0;
3722
3723                              bool addToSwapped = false;
3724
3725                              if (alpha < 0.0) {
3726                                   //at upper bound
3727                                   if (value >= 0.0) {
3728                                        addToSwapped = true;
3729#if TRYBIAS==1
3730                                        badDj = -CoinMax(dj_[iSequence], 0.0);
3731#elif TRYBIAS==2
3732                                        badDj = -dj_[iSequence];
3733#else
3734                                        badDj = -dj_[iSequence] - dualTolerance_;
3735#endif
3736                                   }
3737                              } else {
3738                                   // at lower bound
3739                                   if (value <= 0.0) {
3740                                        addToSwapped = true;
3741#if TRYBIAS==1
3742                                        badDj = CoinMin(dj_[iSequence], 0.0);
3743#elif TRYBIAS==2
3744                                        badDj = dj_[iSequence];
3745#else
3746                                        badDj = dj_[iSequence] - dualTolerance_;
3747#endif
3748                                   }
3749                              }
3750                              if (!addToSwapped) {
3751                                   // add to list of remaining
3752                                   spare2[numberRemaining] = alpha;
3753                                   index2[numberRemaining++] = iSequence;
3754                              } else {
3755                                   // add to list of swapped
3756                                   spare2[--numberPossiblySwapped] = alpha;
3757                                   index2[numberPossiblySwapped] = iSequence;
3758                                   // select if largest pivot
3759                                   bool take = false;
3760                                   double absAlpha = fabs(alpha);
3761#ifdef DUBIOUS_WEIGHTS
3762                                   // User could do anything to break ties here
3763                                   double weight;
3764                                   if (dubiousWeights)
3765                                        weight = dubiousWeights[iSequence];
3766                                   else
3767                                        weight = 1.0;
3768                                   weight += randomNumberGenerator_.randomDouble() * 1.0e-2;
3769                                   if (absAlpha > 2.0 * bestPivot) {
3770                                        take = true;
3771                                   } else if (absAlpha > largestPivot) {
3772                                        // could multiply absAlpha and weight
3773                                        if (weight * bestPivot < bestWeight * absAlpha)
3774                                             take = true;
3775                                   }
3776#else
3777                                   if (absAlpha > bestPivot)
3778                                        take = true;
3779#endif
3780#ifdef MORE_CAREFUL
3781                                   if (absAlpha < acceptablePivot && upperTheta < 1.0e20) {
3782                                        if (alpha < 0.0) {
3783                                             //at upper bound
3784                                             if (value > dualTolerance_) {
3785                                                  double gap = upper_[iSequence] - lower_[iSequence];
3786                                                  if (gap < 1.0e20)
3787                                                       sumBadPivots += value * gap;
3788                                                  else
3789                                                       sumBadPivots += 1.0e20;
3790                                                  //printf("bad %d alpha %g dj at upper %g\n",
3791                                                  //     iSequence,alpha,value);
3792                                             }
3793                                        } else {
3794                                             //at lower bound
3795                                             if (value < -dualTolerance_) {
3796                                                  double gap = upper_[iSequence] - lower_[iSequence];
3797                                                  if (gap < 1.0e20)
3798                                                       sumBadPivots -= value * gap;
3799                                                  else
3800                                                       sumBadPivots += 1.0e20;
3801                                                  //printf("bad %d alpha %g dj at lower %g\n",
3802                                                  //     iSequence,alpha,value);
3803                                             }
3804                                        }
3805                                   }
3806#endif
3807#ifdef FORCE_FOLLOW
3808                                   if (iSequence == force_in) {
3809                                        printf("taking %d - alpha %g best %g\n", force_in, absAlpha, largestPivot);
3810                                        take = true;
3811                                   }
3812#endif
3813                                   if (take) {
3814                                        sequenceIn_ = numberPossiblySwapped;
3815                                        bestPivot =  absAlpha;
3816                                        theta_ = dj_[iSequence] / alpha;
3817                                        largestPivot = CoinMax(largestPivot, 0.5 * bestPivot);
3818#ifdef DUBIOUS_WEIGHTS
3819                                        bestWeight = weight;
3820#endif
3821                                        //printf(" taken seq %d alpha %g weight %d\n",
3822                                        //   iSequence,absAlpha,dubiousWeights[iSequence]);
3823                                   } else {
3824                                        //printf(" not taken seq %d alpha %g weight %d\n",
3825                                        //   iSequence,absAlpha,dubiousWeights[iSequence]);
3826                                   }
3827                                   double range = upper_[iSequence] - lower_[iSequence];
3828                                   thruThis += range * fabs(alpha);
3829                                   increaseInThis += badDj * range;
3830                              }
3831                         }
3832                         marker[1-iFlip][0] = CoinMax(marker[1-iFlip][0], numberRemaining);
3833                         marker[1-iFlip][1] = CoinMin(marker[1-iFlip][1], numberPossiblySwapped);
3834#ifdef MORE_CAREFUL
3835                         // If we have done pivots and things look bad set alpha_ 0.0 to force factorization
3836                         if (sumBadPivots > 1.0e4) {
3837                              if (handler_->logLevel() > 1)
3838                                   *handler_ << "maybe forcing re-factorization - sum " << sumBadPivots << " " << factorization_->pivots() << " pivots" << CoinMessageEol;
3839                              if(factorization_->pivots() > 3) {
3840                                   badSumPivots = true;
3841                                   break;
3842                              }
3843                         }
3844#endif
3845                         swapped[1-iFlip] = numberPossiblySwapped;
3846                         interesting[1-iFlip] = numberRemaining;
3847                         // If we stop now this will be increase in objective (I think)
3848                         double increase = (fabs(dualOut_) - totalThru) * theta_;
3849                         increase += increaseInObjective;
3850                         if (theta_ < 0.0)
3851                              thruThis += fabs(dualOut_); // force using this one
3852                         if (increaseInObjective < 0.0 && increase < 0.0 && lastSequence >= 0) {
3853                              // back
3854                              // We may need to be more careful - we could do by
3855                              // switch so we always do fine grained?
3856                              bestPivot = 0.0;
3857                         } else {
3858                              // add in
3859                              totalThru += thruThis;
3860                              increaseInObjective += increaseInThis;
3861                         }
3862                         if (bestPivot < 0.1 * bestEverPivot &&
3863                                   bestEverPivot > 1.0e-6 &&
3864                                   (bestPivot < 1.0e-3 || totalThru * 2.0 > fabs(dualOut_))) {
3865                              // back to previous one
3866                              sequenceIn_ = lastSequence;
3867                              // swap regions
3868                              iFlip = 1 - iFlip;
3869                              break;
3870                         } else if (sequenceIn_ == -1 && upperTheta > largeValue_) {
3871                              if (lastPivot > acceptablePivot) {
3872                                   // back to previous one
3873                                   sequenceIn_ = lastSequence;
3874                                   // swap regions
3875                                   iFlip = 1 - iFlip;
3876                              } else {
3877                                   // can only get here if all pivots too small
3878                              }
3879                              break;
3880                         } else if (totalThru >= fabs(dualOut_)) {
3881                              modifyCosts = true; // fine grain - we can modify costs
3882                              break; // no point trying another loop
3883                         } else {
3884                              lastSequence = sequenceIn_;
3885                              if (bestPivot > bestEverPivot)
3886                                   bestEverPivot = bestPivot;
3887                              iFlip = 1 - iFlip;
3888                              modifyCosts = true; // fine grain - we can modify costs
3889                         }
3890                    }
3891                    if (iTry == MAXTRY)
3892                         iFlip = 1 - iFlip; // flip back
3893                    break;
3894               } else {
3895                    // skip this lot
3896                    if (bestPivot > 1.0e-3 || bestPivot > bestEverPivot) {
3897                         bestEverPivot = bestPivot;
3898                         lastSequence = bestSequence;
3899                    } else {
3900                         // keep old swapped
3901                         CoinMemcpyN(array[iFlip] + swapped[iFlip],
3902                                     numberColumns_ - swapped[iFlip], array[1-iFlip] + swapped[iFlip]);
3903                         CoinMemcpyN(indices[iFlip] + swapped[iFlip],
3904                                     numberColumns_ - swapped[iFlip], indices[1-iFlip] + swapped[iFlip]);
3905                         marker[1-iFlip][1] = CoinMin(marker[1-iFlip][1], swapped[iFlip]);
3906                         swapped[1-iFlip] = swapped[iFlip];
3907                    }
3908                    increaseInObjective += increaseInThis;
3909                    iFlip = 1 - iFlip; // swap regions
3910                    tentativeTheta = 2.0 * upperTheta;
3911                    totalThru += thruThis;
3912               }
3913          }
3914
3915          // can get here without sequenceIn_ set but with lastSequence
3916          if (sequenceIn_ < 0 && lastSequence >= 0) {
3917               // back to previous one
3918               sequenceIn_ = lastSequence;
3919               // swap regions
3920               iFlip = 1 - iFlip;
3921          }
3922
3923#define MINIMUMTHETA 1.0e-18
3924          // Movement should be minimum for anti-degeneracy - unless
3925          // fixed variable out
3926          double minimumTheta;
3927          if (upperOut_ > lowerOut_)
3928               minimumTheta = MINIMUMTHETA;
3929          else
3930               minimumTheta = 0.0;
3931          if (sequenceIn_ >= 0) {
3932               // at this stage sequenceIn_ is just pointer into index array
3933               // flip just so we can use iFlip
3934               iFlip = 1 - iFlip;
3935               spare = array[iFlip];
3936               index = indices[iFlip];
3937               double oldValue;
3938               alpha_ = spare[sequenceIn_];
3939               sequenceIn_ = indices[iFlip][sequenceIn_];
3940               oldValue = dj_[sequenceIn_];
3941               theta_ = CoinMax(oldValue / alpha_, 0.0);
3942               if (theta_ < minimumTheta && fabs(alpha_) < 1.0e5 && 1) {
3943                    // can't pivot to zero
3944#if 0
3945                    if (oldValue - minimumTheta*alpha_ >= -dualTolerance_) {
3946                         theta_ = minimumTheta;
3947                    } else if (oldValue - minimumTheta*alpha_ >= -newTolerance) {
3948                         theta_ = minimumTheta;
3949                         thisIncrease = true;
3950                    } else {
3951                         theta_ = CoinMax((oldValue + newTolerance) / alpha_, 0.0);
3952                         thisIncrease = true;
3953                    }
3954#else
3955                    theta_ = minimumTheta;
3956#endif
3957               }
3958               // may need to adjust costs so all dual feasible AND pivoted is exactly 0
3959               //int costOffset = numberRows_+numberColumns_;
3960               if (modifyCosts && !badSumPivots) {
3961                    int i;
3962                    for (i = numberColumns_ - 1; i >= swapped[iFlip]; i--) {
3963                         int iSequence = index[i];
3964                         double alpha = spare[i];
3965                         double value = dj_[iSequence] - theta_ * alpha;
3966
3967                         // can't be free here
3968
3969                         if (alpha < 0.0) {
3970                              //at upper bound
3971                              if (value > dualTolerance_) {
3972                                   //thisIncrease = true;
3973#if CLP_CAN_HAVE_ZERO_OBJ<2
3974#define MODIFYCOST 2
3975#endif
3976#if MODIFYCOST
3977                                   // modify cost to hit new tolerance
3978                                   double modification = alpha * theta_ - dj_[iSequence]
3979                                                         + newTolerance;
3980                                   if ((specialOptions_&(2048 + 4096 + 16384)) != 0) {
3981                                        if ((specialOptions_ & 16384) != 0) {
3982                                             if (fabs(modification) < 1.0e-8)
3983                                                  modification = 0.0;
3984                                        } else if ((specialOptions_ & 2048) != 0) {
3985                                             if (fabs(modification) < 1.0e-10)
3986                                                  modification = 0.0;
3987                                        } else {
3988                                             if (fabs(modification) < 1.0e-12)
3989                                                  modification = 0.0;
3990                                        }
3991                                   }
3992                                   dj_[iSequence] += modification;
3993                                   cost_[iSequence] +=  modification;
3994                                   if (modification)
3995                                        numberChanged_ ++; // Say changed costs
3996                                   //cost_[iSequence+costOffset] += modification; // save change
3997#endif
3998                              }
3999                         } else {
4000                              // at lower bound
4001                              if (-value > dualTolerance_) {
4002                                   //thisIncrease = true;
4003#if MODIFYCOST
4004                                   // modify cost to hit new tolerance
4005                                   double modification = alpha * theta_ - dj_[iSequence]
4006                                                         - newTolerance;
4007                                   //modification = CoinMax(modification,-dualTolerance_);
4008                                   //assert (fabs(modification)<1.0e-7);
4009                                   if ((specialOptions_&(2048 + 4096)) != 0) {
4010                                        if ((specialOptions_ & 2048) != 0) {
4011                                             if (fabs(modification) < 1.0e-10)
4012                                                  modification = 0.0;
4013                                        } else {
4014                                             if (fabs(modification) < 1.0e-12)
4015                                                  modification = 0.0;
4016                                        }
4017                                   }
4018                                   dj_[iSequence] += modification;
4019                                   cost_[iSequence] +=  modification;
4020                                   if (modification)
4021                                        numberChanged_ ++; // Say changed costs
4022                                   //cost_[iSequence+costOffset] += modification; // save change
4023#endif
4024                              }
4025                         }
4026                    }
4027               }
4028          }
4029     }
4030
4031#ifdef MORE_CAREFUL
4032     // If we have done pivots and things look bad set alpha_ 0.0 to force factorization
4033     if ((badSumPivots ||
4034               fabs(theta_ * badFree) > 10.0 * dualTolerance_) && factorization_->pivots()) {
4035          if (handler_->logLevel() > 1)
4036               *handler_ << "forcing re-factorization" << CoinMessageEol;
4037          //printf("badSumPivots %g theta_ %g badFree %g\n",badSumPivots,theta_,badFree);
4038          sequenceIn_ = -1;
4039          acceptablePivot_=-acceptablePivot_;
4040     }
4041#endif
4042     if (sequenceIn_ >= 0) {
4043          lowerIn_ = lower_[sequenceIn_];
4044          upperIn_ = upper_[sequenceIn_];
4045          valueIn_ = solution_[sequenceIn_];
4046          dualIn_ = dj_[sequenceIn_];
4047
4048          if (numberTimesOptimal_) {
4049               // can we adjust cost back closer to original
4050               //*** add coding
4051          }
4052#if MODIFYCOST>1
4053          // modify cost to hit zero exactly
4054          // so (dualIn_+modification)==theta_*alpha_
4055          double modification = theta_ * alpha_ - dualIn_;
4056          // But should not move objective too much ??
4057#define DONT_MOVE_OBJECTIVE
4058#ifdef DONT_MOVE_OBJECTIVE
4059          double moveObjective = fabs(modification * solution_[sequenceIn_]);
4060          double smallMove = CoinMax(fabs(objectiveValue_), 1.0e-3);
4061          if (moveObjective > smallMove) {
4062               if (handler_->logLevel() > 1)
4063                    printf("would move objective by %g - original mod %g sol value %g\n", moveObjective,
4064                           modification, solution_[sequenceIn_]);
4065               modification *= smallMove / moveObjective;
4066          }
4067#endif
4068          if (badSumPivots)
4069               modification = 0.0;
4070          if ((specialOptions_&(2048 + 4096)) != 0) {
4071               if ((specialOptions_ & 16384) != 0) {
4072                    // in fast dual
4073                    if (fabs(modification) < 1.0e-7)
4074                         modification = 0.0;
4075               } else if ((specialOptions_ & 2048) != 0) {
4076                    if (fabs(modification) < 1.0e-10)
4077                         modification = 0.0;
4078               } else {
4079                    if (fabs(modification) < 1.0e-12)
4080                         modification = 0.0;
4081               }
4082          }
4083          dualIn_ += modification;
4084          dj_[sequenceIn_] = dualIn_;
4085          cost_[sequenceIn_] += modification;
4086          if (modification)
4087               numberChanged_ ++; // Say changed costs
4088          //int costOffset = numberRows_+numberColumns_;
4089          //cost_[sequenceIn_+costOffset] += modification; // save change
4090          //assert (fabs(modification)<1.0e-6);
4091#ifdef CLP_DEBUG
4092          if ((handler_->logLevel() & 32) && fabs(modification) > 1.0e-15)
4093               printf("exact %d new cost %g, change %g\n", sequenceIn_,
4094                      cost_[sequenceIn_], modification);
4095#endif
4096#endif
4097
4098          if (alpha_ < 0.0) {
4099               // as if from upper bound
4100               directionIn_ = -1;
4101               upperIn_ = valueIn_;
4102          } else {
4103               // as if from lower bound
4104               directionIn_ = 1;
4105               lowerIn_ = valueIn_;
4106          }
4107     } else {
4108          // no pivot
4109          bestPossible = 0.0;
4110          alpha_ = 0.0;
4111     }
4112     //if (thisIncrease)
4113     //dualTolerance_+= 1.0e-6*dblParam_[ClpDualTolerance];
4114
4115     // clear arrays
4116
4117     for (i = 0; i < 2; i++) {
4118          CoinZeroN(array[i], marker[i][0]);
4119          CoinZeroN(array[i] + marker[i][1], numberColumns_ - marker[i][1]);
4120     }
4121     return bestPossible;
4122}
4123#ifdef CLP_ALL_ONE_FILE
4124#undef MAXTRY
4125#endif
4126/* Checks if tentative optimal actually means unbounded
4127   Returns -3 if not, 2 if is unbounded */
4128int
4129ClpSimplexDual::checkUnbounded(CoinIndexedVector * ray,
4130                               CoinIndexedVector * spare,
4131                               double changeCost)
4132{
4133     int status = 2; // say unbounded
4134     factorization_->updateColumn(spare, ray);
4135     // get reduced cost
4136     int i;
4137     int number = ray->getNumElements();
4138     int * index = ray->getIndices();
4139     double * array = ray->denseVector();
4140     for (i = 0; i < number; i++) {
4141          int iRow = index[i];
4142          int iPivot = pivotVariable_[iRow];
4143          changeCost -= cost(iPivot) * array[iRow];
4144     }
4145     double way;
4146     if (changeCost > 0.0) {
4147          //try going down
4148          way = 1.0;
4149     } else if (changeCost < 0.0) {
4150          //try going up
4151          way = -1.0;
4152     } else {
4153#ifdef CLP_DEBUG
4154          printf("can't decide on up or down\n");
4155#endif
4156          way = 0.0;
4157          status = -3;
4158     }
4159     double movement = 1.0e10 * way; // some largish number
4160     double zeroTolerance = 1.0e-14 * dualBound_;
4161     for (i = 0; i < number; i++) {
4162          int iRow = index[i];
4163          int iPivot = pivotVariable_[iRow];
4164          double arrayValue = array[iRow];
4165          if (fabs(arrayValue) < zeroTolerance)
4166               arrayValue = 0.0;
4167          double newValue = solution(iPivot) + movement * arrayValue;
4168          if (newValue > upper(iPivot) + primalTolerance_ ||
4169                    newValue < lower(iPivot) - primalTolerance_)
4170               status = -3; // not unbounded
4171     }
4172     if (status == 2) {
4173          // create ray
4174          delete [] ray_;
4175          ray_ = new double [numberColumns_];
4176          CoinZeroN(ray_, numberColumns_);
4177          for (i = 0; i < number; i++) {
4178               int iRow = index[i];
4179               int iPivot = pivotVariable_[iRow];
4180               double arrayValue = array[iRow];
4181               if (iPivot < numberColumns_ && fabs(arrayValue) >= zeroTolerance)
4182                    ray_[iPivot] = way * array[iRow];
4183          }
4184     }
4185     ray->clear();
4186     return status;
4187}
4188//static int count_status=0;
4189//static double obj_status=0.0;
4190//static int check_status=123456789;//41754;
4191//static int count_alpha=0;
4192/* Checks if finished.  Updates status */
4193void
4194ClpSimplexDual::statusOfProblemInDual(int & lastCleaned, int type,
4195                                      double * givenDuals, ClpDataSave & saveData,
4196                                      int ifValuesPass)
4197{
4198#ifdef CLP_INVESTIGATE_SERIAL
4199     if (z_thinks > 0 && z_thinks < 2)
4200          z_thinks += 2;
4201#endif
4202     bool arraysNotCreated = (type==0);
4203     // If lots of iterations then adjust costs if large ones
4204     if (numberIterations_ > 4 * (numberRows_ + numberColumns_) && objectiveScale_ == 1.0) {
4205          double largest = 0.0;
4206          for (int i = 0; i < numberRows_; i++) {
4207               int iColumn = pivotVariable_[i];
4208               largest = CoinMax(largest, fabs(cost_[iColumn]));
4209          }
4210          if (largest > 1.0e6) {
4211               objectiveScale_ = 1.0e6 / largest;
4212               for (int i = 0; i < numberRows_ + numberColumns_; i++)
4213                    cost_[i] *= objectiveScale_;
4214          }
4215     }
4216     int numberPivots = factorization_->pivots();
4217     double realDualInfeasibilities = 0.0;
4218     if (type == 2) {
4219          if (alphaAccuracy_ != -1.0)
4220               alphaAccuracy_ = -2.0;
4221          // trouble - restore solution
4222          CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4223          CoinMemcpyN(savedSolution_ + numberColumns_ ,
4224                      numberRows_, rowActivityWork_);
4225          CoinMemcpyN(savedSolution_ ,
4226                      numberColumns_, columnActivityWork_);
4227          // restore extra stuff
4228          int dummy;
4229          matrix_->generalExpanded(this, 6, dummy);
4230          forceFactorization_ = 1; // a bit drastic but ..
4231          changeMade_++; // say something changed
4232          // get correct bounds on all variables
4233          resetFakeBounds(0);
4234     }
4235     int tentativeStatus = problemStatus_;
4236     double changeCost;
4237     bool unflagVariables = true;
4238     bool weightsSaved = false;
4239     bool weightsSaved2 = numberIterations_ && !numberPrimalInfeasibilities_;
4240     int dontFactorizePivots = dontFactorizePivots_;
4241     if (type == 3) {
4242          type = 1;
4243          dontFactorizePivots = 1;
4244     }
4245     if (alphaAccuracy_ < 0.0 || !numberPivots || alphaAccuracy_ > 1.0e4 || numberPivots > 20) {
4246          if (problemStatus_ > -3 || numberPivots > dontFactorizePivots) {
4247               // factorize
4248               // later on we will need to recover from singularities
4249               // also we could skip if first time
4250               // save dual weights
4251               dualRowPivot_->saveWeights(this, 1);
4252               weightsSaved = true;
4253               if (type) {
4254                    // is factorization okay?
4255                    if (internalFactorize(1)) {
4256                         // no - restore previous basis
4257                         unflagVariables = false;
4258                         assert (type == 1);
4259                         changeMade_++; // say something changed
4260                         // Keep any flagged variables
4261                         int i;
4262                         for (i = 0; i < numberRows_ + numberColumns_; i++) {
4263                              if (flagged(i))
4264                                   saveStatus_[i] |= 64; //say flagged
4265                         }
4266                         CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4267                         CoinMemcpyN(savedSolution_ + numberColumns_ ,
4268                                     numberRows_, rowActivityWork_);
4269                         CoinMemcpyN(savedSolution_ ,
4270                                     numberColumns_, columnActivityWork_);
4271                         // restore extra stuff
4272                         int dummy;
4273                         matrix_->generalExpanded(this, 6, dummy);
4274                         // get correct bounds on all variables
4275                         resetFakeBounds(1);
4276                         // need to reject something
4277                         char x = isColumn(sequenceOut_) ? 'C' : 'R';
4278                         handler_->message(CLP_SIMPLEX_FLAG, messages_)
4279                                   << x << sequenceWithin(sequenceOut_)
4280                                   << CoinMessageEol;
4281#ifdef COIN_DEVELOP
4282                         printf("flag d\n");
4283#endif
4284                         setFlagged(sequenceOut_);
4285                         progress_.clearBadTimes();
4286
4287                         // Go to safe
4288                         factorization_->pivotTolerance(0.99);
4289                         forceFactorization_ = 1; // a bit drastic but ..
4290                         type = 2;
4291                         //assert (internalFactorize(1)==0);
4292                         if (internalFactorize(1)) {
4293                              CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4294                              CoinMemcpyN(savedSolution_ + numberColumns_ ,
4295                                          numberRows_, rowActivityWork_);
4296                              CoinMemcpyN(savedSolution_ ,
4297                                          numberColumns_, columnActivityWork_);
4298                              // restore extra stuff
4299                              int dummy;
4300                              matrix_->generalExpanded(this, 6, dummy);
4301                              // debug
4302                              int returnCode = internalFactorize(1);
4303                              while (returnCode) {
4304                                   // ouch
4305                                   // switch off dense
4306                                   int saveDense = factorization_->denseThreshold();
4307                                   factorization_->setDenseThreshold(0);
4308                                   // Go to safe
4309                                   factorization_->pivotTolerance(0.99);
4310                                   // make sure will do safe factorization
4311                                   pivotVariable_[0] = -1;
4312                                   returnCode = internalFactorize(2);
4313                                   factorization_->setDenseThreshold(saveDense);
4314                              }
4315                              // get correct bounds on all variables
4316                              resetFakeBounds(1);
4317                         }
4318                    }
4319               }
4320               if (problemStatus_ != -4 || numberPivots > 10)
4321                    problemStatus_ = -3;
4322          }
4323     } else {
4324          //printf("testing with accuracy of %g and status of %d\n",alphaAccuracy_,problemStatus_);
4325          //count_alpha++;
4326          //if ((count_alpha%5000)==0)
4327          //printf("count alpha %d\n",count_alpha);
4328     }
4329     if(progress_.infeasibility_[0]<1.0e-1 &&
4330        primalTolerance_==1.0e-7&&progress_.iterationNumber_[0]>0&&
4331        progress_.iterationNumber_[CLP_PROGRESS-1]-progress_.iterationNumber_[0]>25) {
4332       // default - so user did not set
4333       int iP;
4334       double minAverage=COIN_DBL_MAX;
4335       double maxAverage=0.0;
4336       for (iP=0;iP<CLP_PROGRESS;iP++) {
4337         int n=progress_.numberInfeasibilities_[iP];
4338         if (!n) {
4339           break;
4340         } else {
4341           double average=progress_.infeasibility_[iP];
4342           if (average>0.1)
4343             break;
4344           average /= static_cast<double>(n);
4345           minAverage=CoinMin(minAverage,average);
4346           maxAverage=CoinMax(maxAverage,average);
4347         }
4348       }
4349       if (iP==CLP_PROGRESS&&minAverage<1.0e-5&&maxAverage<1.0e-3) {
4350         // change tolerance
4351#if CBC_USEFUL_PRINTING>0
4352         printf("CCchanging tolerance\n");
4353#endif
4354         primalTolerance_=1.0e-6;
4355         minimumPrimalTolerance_=primalTolerance_;
4356         dblParam_[ClpPrimalTolerance]=1.0e-6;
4357         moreSpecialOptions_ |= 4194304;
4358       }
4359     }
4360     // at this stage status is -3 or -4 if looks infeasible
4361     // get primal and dual solutions
4362#if 0
4363     {
4364          int numberTotal = numberRows_ + numberColumns_;
4365          double * saveSol = CoinCopyOfArray(solution_, numberTotal);
4366          double * saveDj = CoinCopyOfArray(dj_, numberTotal);
4367          double tolerance = type ? 1.0e-4 : 1.0e-8;
4368          // always if values pass
4369          double saveObj = objectiveValue_;
4370          double sumPrimal = sumPrimalInfeasibilities_;
4371          int numberPrimal = numberPrimalInfeasibilities_;
4372          double sumDual = sumDualInfeasibilities_;
4373          int numberDual = numberDualInfeasibilities_;
4374          gutsOfSolution(givenDuals, NULL);
4375          int j;
4376          double largestPrimal = tolerance;
4377          int iPrimal = -1;
4378          for (j = 0; j < numberTotal; j++) {
4379               double difference = solution_[j] - saveSol[j];
4380               if (fabs(difference) > largestPrimal) {
4381                    iPrimal = j;
4382                    largestPrimal = fabs(difference);
4383               }
4384          }
4385          double largestDual = tolerance;
4386          int iDual = -1;
4387          for (j = 0; j < numberTotal; j++) {
4388               double difference = dj_[j] - saveDj[j];
4389               if (fabs(difference) > largestDual && upper_[j] > lower_[j]) {
4390                    iDual = j;
4391                    largestDual = fabs(difference);
4392               }
4393          }
4394          if (!type) {
4395               if (fabs(saveObj - objectiveValue_) > 1.0e-5 ||
4396                         numberPrimal != numberPrimalInfeasibilities_ || numberPrimal != 1 ||
4397                         fabs(sumPrimal - sumPrimalInfeasibilities_) > 1.0e-5 || iPrimal >= 0 ||
4398                         numberDual != numberDualInfeasibilities_ || numberDual != 0 ||
4399                         fabs(sumDual - sumDualInfeasibilities_) > 1.0e-5 || iDual >= 0)
4400                    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",
4401                           type, numberIterations_, numberPivots,
4402                           numberPrimal, numberPrimalInfeasibilities_, sumPrimal, sumPrimalInfeasibilities_,
4403                           largestPrimal, iPrimal,
4404                           numberDual, numberDualInfeasibilities_, sumDual, sumDualInfeasibilities_,
4405                           largestDual, iDual,
4406                           saveObj, objectiveValue_);
4407          } else {
4408               if (fabs(saveObj - objectiveValue_) > 1.0e-5 ||
4409                         numberPrimalInfeasibilities_ || iPrimal >= 0 ||
4410                         numberDualInfeasibilities_ || iDual >= 0)
4411                    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",
4412                           type, numberIterations_, numberPivots,
4413                           numberPrimal, numberPrimalInfeasibilities_, sumPrimal, sumPrimalInfeasibilities_,
4414                           largestPrimal, iPrimal,
4415                           numberDual, numberDualInfeasibilities_, sumDual, sumDualInfeasibilities_,
4416                           largestDual, iDual,
4417                           saveObj, objectiveValue_);
4418          }
4419          delete [] saveSol;
4420          delete [] saveDj;
4421     }
4422#else
4423     if (type || ifValuesPass)
4424          gutsOfSolution(givenDuals, NULL);
4425#endif
4426     // If bad accuracy treat as singular
4427     if ((largestPrimalError_ > 1.0e15 || largestDualError_ > 1.0e15) && numberIterations_) {
4428          // restore previous basis
4429          unflagVariables = false;
4430          changeMade_++; // say something changed
4431          // Keep any flagged variables
4432          int i;
4433          for (i = 0; i < numberRows_ + numberColumns_; i++) {
4434               if (flagged(i))
4435                    saveStatus_[i] |= 64; //say flagged
4436          }
4437          CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4438          CoinMemcpyN(savedSolution_ + numberColumns_ ,
4439                      numberRows_, rowActivityWork_);
4440          CoinMemcpyN(savedSolution_ ,
4441                      numberColumns_, columnActivityWork_);
4442          // restore extra stuff
4443          int dummy;
4444          matrix_->generalExpanded(this, 6, dummy);
4445          // get correct bounds on all variables
4446          resetFakeBounds(1);
4447          // need to reject something
4448          char x = isColumn(sequenceOut_) ? 'C' : 'R';
4449          handler_->message(CLP_SIMPLEX_FLAG, messages_)
4450                    << x << sequenceWithin(sequenceOut_)
4451                    << CoinMessageEol;
4452#ifdef COIN_DEVELOP
4453          printf("flag e\n");
4454#endif
4455          setFlagged(sequenceOut_);
4456          progress_.clearBadTimes();
4457
4458          // Go to safer
4459          double newTolerance = CoinMin(1.1 * factorization_->pivotTolerance(), 0.99);
4460          factorization_->pivotTolerance(newTolerance);
4461          forceFactorization_ = 1; // a bit drastic but ..
4462          if (alphaAccuracy_ != -1.0)
4463               alphaAccuracy_ = -2.0;
4464          type = 2;
4465          //assert (internalFactorize(1)==0);
4466          if (internalFactorize(1)) {
4467               CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4468               CoinMemcpyN(savedSolution_ + numberColumns_ ,
4469                           numberRows_, rowActivityWork_);
4470               CoinMemcpyN(savedSolution_ ,
4471                           numberColumns_, columnActivityWork_);
4472               // restore extra stuff
4473               int dummy;
4474               matrix_->generalExpanded(this, 6, dummy);
4475               // debug
4476               int returnCode = internalFactorize(1);
4477               while (returnCode) {
4478                    // ouch
4479                    // switch off dense
4480                    int saveDense = factorization_->denseThreshold();
4481                    factorization_->setDenseThreshold(0);
4482                    // Go to safe
4483                    factorization_->pivotTolerance(0.99);
4484                    // make sure will do safe factorization
4485                    pivotVariable_[0] = -1;
4486                    returnCode = internalFactorize(2);
4487                    factorization_->setDenseThreshold(saveDense);
4488               }
4489               // get correct bounds on all variables
4490               resetFakeBounds(1);
4491          }
4492          // get primal and dual solutions
4493          gutsOfSolution(givenDuals, NULL);
4494     } else if (goodAccuracy()) {
4495          // Can reduce tolerance
4496          double newTolerance = CoinMax(0.995 * factorization_->pivotTolerance(), saveData.pivotTolerance_);
4497          factorization_->pivotTolerance(newTolerance);
4498     }
4499     bestObjectiveValue_ = CoinMax(bestObjectiveValue_,
4500                                   objectiveValue_ - bestPossibleImprovement_);
4501     bool reallyBadProblems = false;
4502     // Double check infeasibility if no action
4503     if (progress_.lastIterationNumber(0) == numberIterations_) {
4504          if (dualRowPivot_->looksOptimal()) {
4505               numberPrimalInfeasibilities_ = 0;
4506               sumPrimalInfeasibilities_ = 0.0;
4507          }
4508#if 1
4509     } else {
4510          double thisObj = objectiveValue_ - bestPossibleImprovement_;
4511#ifdef CLP_INVESTIGATE
4512          assert (bestPossibleImprovement_ > -1000.0 && objectiveValue_ > -1.0e100);
4513          if (bestPossibleImprovement_)
4514               printf("obj %g add in %g -> %g\n", objectiveValue_, bestPossibleImprovement_,
4515                      thisObj);
4516#endif
4517          double lastObj = progress_.lastObjective(0);
4518#ifndef NDEBUG
4519#ifdef COIN_DEVELOP
4520          resetFakeBounds(-1);
4521#endif
4522#endif
4523#ifdef CLP_REPORT_PROGRESS
4524          ixxxxxx++;
4525          if (ixxxxxx >= ixxyyyy - 4 && ixxxxxx <= ixxyyyy) {
4526               char temp[20];
4527               sprintf(temp, "sol%d.out", ixxxxxx);
4528               printf("sol%d.out\n", ixxxxxx);
4529               FILE * fp = fopen(temp, "w");
4530               int nTotal = numberRows_ + numberColumns_;
4531               for (int i = 0; i < nTotal; i++)
4532                    fprintf(fp, "%d %d %g %g %g %g %g\n",
4533                            i, status_[i], lower_[i], solution_[i], upper_[i], cost_[i], dj_[i]);
4534               fclose(fp);
4535          }
4536#endif
4537          if(!ifValuesPass && firstFree_ < 0) {
4538               double testTol = 5.0e-3;
4539               if (progress_.timesFlagged() > 10) {
4540                    testTol *= pow(2.0, progress_.timesFlagged() - 8);
4541               } else if (progress_.timesFlagged() > 5) {
4542                    testTol *= 5.0;
4543               }
4544               if (lastObj > thisObj +
4545                         testTol*(fabs(thisObj) + fabs(lastObj)) + testTol) {
4546                    int maxFactor = factorization_->maximumPivots();
4547                    if ((specialOptions_ & 1048576) == 0) {
4548                         if (progress_.timesFlagged() > 10)
4549                              progress_.incrementReallyBadTimes();
4550                         if (maxFactor > 10 - 9) {
4551#ifdef COIN_DEVELOP
4552                              printf("lastobj %g thisobj %g\n", lastObj, thisObj);
4553#endif
4554                              //if (forceFactorization_<0)
4555                              //forceFactorization_= maxFactor;
4556                              //forceFactorization_ = CoinMax(1,(forceFactorization_>>1));
4557                              if ((progressFlag_ & 4) == 0 && lastObj < thisObj + 1.0e4 &&
4558                                        largestPrimalError_ < 1.0e2) {
4559                                   // Just save costs
4560                                   // save extra copy of cost_
4561                                   int nTotal = numberRows_ + numberColumns_;
4562                                   double * temp = new double [2*nTotal];
4563                                   memcpy(temp, cost_, nTotal * sizeof(double));
4564                                   memcpy(temp + nTotal, cost_, nTotal * sizeof(double));
4565                                   delete [] cost_;
4566                                   cost_ = temp;
4567                                   objectiveWork_ = cost_;
4568                                   rowObjectiveWork_ = cost_ + numberColumns_;
4569                                   progressFlag_ |= 4;
4570                              } else {
4571                                   forceFactorization_ = 1;
4572#ifdef COIN_DEVELOP
4573                                   printf("Reducing factorization frequency - bad backwards\n");
4574#endif
4575#if 1
4576                                   unflagVariables = false;
4577                                   changeMade_++; // say something changed
4578                                   int nTotal = numberRows_ + numberColumns_;
4579                                   CoinMemcpyN(saveStatus_, nTotal, status_);
4580                                   CoinMemcpyN(savedSolution_ + numberColumns_ ,
4581                                               numberRows_, rowActivityWork_);
4582                                   CoinMemcpyN(savedSolution_ ,
4583                                               numberColumns_, columnActivityWork_);
4584                                   if ((progressFlag_ & 4) == 0) {
4585                                        // save extra copy of cost_
4586                                        double * temp = new double [2*nTotal];
4587                                        memcpy(temp, cost_, nTotal * sizeof(double));
4588                                        memcpy(temp + nTotal, cost_, nTotal * sizeof(double));
4589                                        delete [] cost_;
4590                                        cost_ = temp;
4591                                        objectiveWork_ = cost_;
4592                                        rowObjectiveWork_ = cost_ + numberColumns_;
4593                                        progressFlag_ |= 4;
4594                                   } else {
4595                                        memcpy(cost_, cost_ + nTotal, nTotal * sizeof(double));
4596                                   }
4597                                   // restore extra stuff
4598                                   int dummy;
4599                                   matrix_->generalExpanded(this, 6, dummy);
4600                                   double pivotTolerance = factorization_->pivotTolerance();
4601                                   if(pivotTolerance < 0.2)
4602                                        factorization_->pivotTolerance(0.2);
4603                                   else if(progress_.timesFlagged() > 2)
4604                                        factorization_->pivotTolerance(CoinMin(pivotTolerance * 1.1, 0.99));
4605                                   if (alphaAccuracy_ != -1.0)
4606                                        alphaAccuracy_ = -2.0;
4607                                   if (internalFactorize(1)) {
4608                                        CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4609                                        CoinMemcpyN(savedSolution_ + numberColumns_ ,
4610                                                    numberRows_, rowActivityWork_);
4611                                        CoinMemcpyN(savedSolution_ ,
4612                                                    numberColumns_, columnActivityWork_);
4613                                        // restore extra stuff
4614                                        int dummy;
4615                                        matrix_->generalExpanded(this, 6, dummy);
4616                                        // debug
4617                                        int returnCode = internalFactorize(1);
4618                                        while (returnCode) {
4619                                             // ouch
4620                                             // switch off dense
4621                                             int saveDense = factorization_->denseThreshold();
4622                                             factorization_->setDenseThreshold(0);
4623                                             // Go to safe
4624                                             factorization_->pivotTolerance(0.99);
4625                                             // make sure will do safe factorization
4626                                             pivotVariable_[0] = -1;
4627                                             returnCode = internalFactorize(2);
4628                                             factorization_->setDenseThreshold(saveDense);
4629                                        }
4630                                   }
4631                                   resetFakeBounds(0);
4632                                   type = 2; // so will restore weights
4633                                   // get primal and dual solutions
4634                                   gutsOfSolution(givenDuals, NULL);
4635                                   if (numberPivots < 2) {
4636                                        // need to reject something
4637                                        char x = isColumn(sequenceOut_) ? 'C' : 'R';
4638                                        handler_->message(CLP_SIMPLEX_FLAG, messages_)
4639                                                  << x << sequenceWithin(sequenceOut_)
4640                                                  << CoinMessageEol;
4641#ifdef COIN_DEVELOP
4642                                        printf("flag d\n");
4643#endif
4644                                        setFlagged(sequenceOut_);
4645                                        progress_.clearBadTimes();
4646                                        progress_.incrementTimesFlagged();
4647                                   }
4648                                   if (numberPivots < 10)
4649                                        reallyBadProblems = true;
4650#ifdef COIN_DEVELOP
4651                                   printf("obj now %g\n", objectiveValue_);
4652#endif
4653                                   progress_.modifyObjective(objectiveValue_
4654                                                             - bestPossibleImprovement_);
4655#endif
4656                              }
4657                         }
4658                    } else {
4659                         // in fast dual give up
4660#ifdef COIN_DEVELOP
4661                         printf("In fast dual?\n");
4662#endif
4663                         problemStatus_ = 3;
4664                    }
4665               } else if (lastObj < thisObj - 1.0e-5 * CoinMax(fabs(thisObj), fabs(lastObj)) - 1.0e-3) {
4666                    numberTimesOptimal_ = 0;
4667               }
4668          }
4669#endif
4670     }
4671     // Up tolerance if looks a bit odd
4672     if (numberIterations_ > CoinMax(1000, numberRows_ >> 4) && (specialOptions_ & 64) != 0) {
4673          if (sumPrimalInfeasibilities_ && sumPrimalInfeasibilities_ < 1.0e5) {
4674               int backIteration = progress_.lastIterationNumber(CLP_PROGRESS - 1);
4675               if (backIteration > 0 && numberIterations_ - backIteration < 9 * CLP_PROGRESS) {
4676                    if (factorization_->pivotTolerance() < 0.9) {
4677                         // up tolerance
4678                         factorization_->pivotTolerance(CoinMin(factorization_->pivotTolerance() * 1.05 + 0.02, 0.91));
4679                         //printf("tol now %g\n",factorization_->pivotTolerance());
4680                         progress_.clearIterationNumbers();
4681                    }
4682               }
4683          }
4684     }
4685     // Check if looping
4686     int loop;
4687     if (!givenDuals && type != 2)
4688          loop = progress_.looping();
4689     else
4690          loop = -1;
4691     if (progress_.reallyBadTimes() > 10) {
4692          problemStatus_ = 10; // instead - try other algorithm
4693#if COIN_DEVELOP>2
4694          printf("returning at %d\n", __LINE__);
4695#endif
4696     }
4697     int situationChanged = 0;
4698     if (loop >= 0) {
4699          problemStatus_ = loop; //exit if in loop
4700          if (!problemStatus_) {
4701               // declaring victory
4702               numberPrimalInfeasibilities_ = 0;
4703               sumPrimalInfeasibilities_ = 0.0;
4704          } else {
4705               problemStatus_ = 10; // instead - try other algorithm
4706#if COIN_DEVELOP>2
4707               printf("returning at %d\n", __LINE__);
4708#endif
4709          }
4710          return;
4711     } else if (loop < -1) {
4712          // something may have changed
4713          gutsOfSolution(NULL, NULL);
4714          situationChanged = 1;
4715     }
4716     // really for free variables in
4717     if((progressFlag_ & 2) != 0) {
4718          situationChanged = 2;
4719     }
4720     progressFlag_ &= (~3); //reset progress flag
4721     if ((progressFlag_ & 4) != 0) {
4722          // save copy of cost_
4723          int nTotal = numberRows_ + numberColumns_;
4724          memcpy(cost_ + nTotal, cost_, nTotal * sizeof(double));
4725     }
4726     /*if (!numberIterations_&&sumDualInfeasibilities_)
4727       printf("OBJ %g sumPinf %g sumDinf %g\n",
4728        objectiveValue(),sumPrimalInfeasibilities_,
4729        sumDualInfeasibilities_);*/
4730     // mark as having gone optimal if looks like it
4731     if (!numberPrimalInfeasibilities_&&
4732         !numberDualInfeasibilities_)
4733       progressFlag_ |= 8;
4734     if (handler_->detail(CLP_SIMPLEX_STATUS, messages_) < 100) {
4735          handler_->message(CLP_SIMPLEX_STATUS, messages_)
4736                    << numberIterations_ << objectiveValue();
4737          handler_->printing(sumPrimalInfeasibilities_ > 0.0)
4738                    << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
4739          handler_->printing(sumDualInfeasibilities_ > 0.0)
4740                    << sumDualInfeasibilities_ << numberDualInfeasibilities_;
4741          handler_->printing(numberDualInfeasibilitiesWithoutFree_
4742                             < numberDualInfeasibilities_)
4743                    << numberDualInfeasibilitiesWithoutFree_;
4744          handler_->message() << CoinMessageEol;
4745     }
4746#if 0
4747     count_status++;
4748     if (!numberIterations_)
4749       obj_status=-1.0e30;
4750     if (objectiveValue()<obj_status-0.01) {
4751       printf("Backward obj at %d from %g to %g\n",
4752              count_status,obj_status,objectiveValue());
4753     }
4754     obj_status=objectiveValue();
4755     if (count_status>=check_status-1) {
4756       printf("Trouble ahead - count_status %d\n",count_status);
4757     }
4758#endif
4759#if 0
4760     printf("IT %d %g %g(%d) %g(%d)\n",
4761            numberIterations_, objectiveValue(),
4762            sumPrimalInfeasibilities_, numberPrimalInfeasibilities_,
4763            sumDualInfeasibilities_, numberDualInfeasibilities_);
4764#endif
4765     double approximateObjective = objectiveValue_;
4766#ifdef CLP_REPORT_PROGRESS
4767     if (ixxxxxx >= ixxyyyy - 4 && ixxxxxx <= ixxyyyy) {
4768          char temp[20];
4769          sprintf(temp, "x_sol%d.out", ixxxxxx);
4770          FILE * fp = fopen(temp, "w");
4771          int nTotal = numberRows_ + numberColumns_;
4772          for (int i = 0; i < nTotal; i++)
4773               fprintf(fp, "%d %d %g %g %g %g %g\n",
4774                       i, status_[i], lower_[i], solution_[i], upper_[i], cost_[i], dj_[i]);
4775          fclose(fp);
4776          if (ixxxxxx == ixxyyyy)
4777               exit(6);
4778     }
4779#endif
4780     realDualInfeasibilities = sumDualInfeasibilities_;
4781     double saveTolerance = dualTolerance_;
4782     // If we need to carry on cleaning variables
4783     if (!numberPrimalInfeasibilities_ && (specialOptions_ & 1024) != 0 && CLEAN_FIXED) {
4784          for (int iRow = 0; iRow < numberRows_; iRow++) {
4785               int iPivot = pivotVariable_[iRow];
4786               if (!flagged(iPivot) && pivoted(iPivot)) {
4787                    // carry on
4788                    numberPrimalInfeasibilities_ = -1;
4789                    sumOfRelaxedPrimalInfeasibilities_ = 1.0;
4790                    sumPrimalInfeasibilities_ = 1.0;
4791                    break;
4792               }
4793          }
4794     }
4795     /* If we are primal feasible and any dual infeasibilities are on
4796        free variables then it is better to go to primal */
4797     if (!numberPrimalInfeasibilities_ && ((!numberDualInfeasibilitiesWithoutFree_ &&
4798                                            numberDualInfeasibilities_)||
4799                                           (moreSpecialOptions_&2097152)!=0))
4800          problemStatus_ = 10;
4801     // dual bound coming in
4802     //double saveDualBound = dualBound_;
4803     bool needCleanFake = false;
4804     while (problemStatus_ <= -3) {
4805          int cleanDuals = 0;
4806          if (situationChanged != 0)
4807               cleanDuals = 1;
4808          int numberChangedBounds = 0;
4809          int doOriginalTolerance = 0;
4810          if ( lastCleaned == numberIterations_)
4811               doOriginalTolerance = 1;
4812          // check optimal
4813          // give code benefit of doubt
4814          if (sumOfRelaxedDualInfeasibilities_ == 0.0 &&
4815                    sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
4816               // say optimal (with these bounds etc)
4817               numberDualInfeasibilities_ = 0;
4818               sumDualInfeasibilities_ = 0.0;
4819               numberPrimalInfeasibilities_ = 0;
4820               sumPrimalInfeasibilities_ = 0.0;
4821          }
4822          //if (dualFeasible()||problemStatus_==-4||(primalFeasible()&&!numberDualInfeasibilitiesWithoutFree_)) {
4823          if (dualFeasible() || problemStatus_ == -4) {
4824               progress_.modifyObjective(objectiveValue_
4825                                         - bestPossibleImprovement_);
4826#ifdef COIN_DEVELOP
4827               if (sumDualInfeasibilities_ || bestPossibleImprovement_)
4828                    printf("improve %g dualinf %g -> %g\n",
4829                           bestPossibleImprovement_, sumDualInfeasibilities_,
4830                           sumDualInfeasibilities_ * dualBound_);
4831#endif
4832               // see if cutoff reached
4833               double limit = 0.0;
4834               getDblParam(ClpDualObjectiveLimit, limit);
4835#if 0
4836               if(fabs(limit) < 1.0e30 && objectiveValue()*optimizationDirection_ >
4837                         limit + 1.0e-7 + 1.0e-8 * fabs(limit) && !numberAtFakeBound()) {
4838                    //looks infeasible on objective
4839                    if (perturbation_ == 101) {
4840                         cleanDuals = 1;
4841                         // Save costs
4842                         int numberTotal = numberRows_ + numberColumns_;
4843                         double * saveCost = CoinCopyOfArray(cost_, numberTotal);
4844                         // make sure fake bounds are back
4845                         changeBounds(1, NULL, changeCost);
4846                         createRim4(false);
4847                         // make sure duals are current
4848                         computeDuals(givenDuals);
4849                         checkDualSolution();
4850                         if(objectiveValue()*optimizationDirection_ >
4851                                   limit + 1.0e-7 + 1.0e-8 * fabs(limit) && !numberDualInfeasibilities_) {
4852                              perturbation_ = 102; // stop any perturbations
4853                              printf("cutoff test succeeded\n");
4854                         } else {
4855                              printf("cutoff test failed\n");
4856                              // put back
4857                              memcpy(cost_, saveCost, numberTotal * sizeof(double));
4858                              // make sure duals are current
4859                              computeDuals(givenDuals);
4860                              checkDualSolution();
4861                              progress_.modifyObjective(-COIN_DBL_MAX);
4862                              problemStatus_ = -1;
4863                         }
4864                         delete [] saveCost;
4865                    }
4866               }
4867#endif
4868               if (primalFeasible() && !givenDuals) {
4869                    // may be optimal - or may be bounds are wrong
4870                    handler_->message(CLP_DUAL_BOUNDS, messages_)
4871                              << dualBound_
4872                              << CoinMessageEol;
4873                    // save solution in case unbounded
4874                    double * saveColumnSolution = NULL;
4875                    double * saveRowSolution = NULL;
4876                    bool inCbc = (specialOptions_ & (0x01000000 | 16384)) != 0;
4877                    if (!inCbc) {
4878                         saveColumnSolution = CoinCopyOfArray(columnActivityWork_, numberColumns_);
4879                         saveRowSolution = CoinCopyOfArray(rowActivityWork_, numberRows_);
4880                    }
4881                    numberChangedBounds = changeBounds(0, rowArray_[3], changeCost);
4882                    if (numberChangedBounds <= 0 && !numberDualInfeasibilities_) {
4883                         //looks optimal - do we need to reset tolerance
4884                         if (perturbation_ == 101) {
4885                              perturbation_ = 102; // stop any perturbations
4886                              cleanDuals = 1;
4887                              // make sure fake bounds are back
4888                              //computeObjectiveValue();
4889                              changeBounds(1, NULL, changeCost);
4890                              //computeObjectiveValue();
4891                              createRim4(false);
4892                              // make sure duals are current
4893                              computeDuals(givenDuals);
4894                              checkDualSolution();
4895                              progress_.modifyObjective(-COIN_DBL_MAX);
4896#define DUAL_TRY_FASTER
4897#ifdef DUAL_TRY_FASTER
4898                              if (numberDualInfeasibilities_) {
4899#endif
4900                                   numberChanged_ = 1; // force something to happen
4901                                   lastCleaned = numberIterations_ - 1;
4902#ifdef DUAL_TRY_FASTER
4903                              } else {
4904                                   //double value = objectiveValue_;
4905                                   computeObjectiveValue(true);
4906                                   //printf("old %g new %g\n",value,objectiveValue_);
4907                                   //numberChanged_=1;
4908                              }
4909#endif
4910                         }
4911                         if (lastCleaned < numberIterations_ && numberTimesOptimal_ < 4 &&
4912                                   (numberChanged_ || (specialOptions_ & 4096) == 0)) {
4913#if CLP_CAN_HAVE_ZERO_OBJ
4914                           if ((specialOptions_&2097152)==0) {
4915#endif
4916                              doOriginalTolerance = 2;
4917                              numberTimesOptimal_++;
4918                              changeMade_++; // say something changed
4919                              if (numberTimesOptimal_ == 1) {
4920                                   dualTolerance_ = dblParam_[ClpDualTolerance];
4921                              } else {
4922                                   if (numberTimesOptimal_ == 2) {
4923                                        // better to have small tolerance even if slower
4924                                        factorization_->zeroTolerance(CoinMin(factorization_->zeroTolerance(), 1.0e-15));
4925                                   }
4926                                   dualTolerance_ = dblParam_[ClpDualTolerance];
4927                                   dualTolerance_ *= pow(2.0, numberTimesOptimal_ - 1);
4928                              }
4929                              cleanDuals = 2; // If nothing changed optimal else primal
4930#if CLP_CAN_HAVE_ZERO_OBJ
4931                           } else {
4932                             // no cost - skip checks
4933                             problemStatus_=0;
4934                           }
4935#endif
4936                         } else {
4937                              problemStatus_ = 0; // optimal
4938                              if (lastCleaned < numberIterations_ && numberChanged_) {
4939                                   handler_->message(CLP_SIMPLEX_GIVINGUP, messages_)
4940                                             << CoinMessageEol;
4941                              }
4942                         }
4943                    } else {
4944                         cleanDuals = 1;
4945                         if (doOriginalTolerance == 1) {
4946                              // check unbounded
4947                              // find a variable with bad dj
4948                              int iSequence;
4949                              int iChosen = -1;
4950                              if (!inCbc) {
4951                                   double largest = 100.0 * primalTolerance_;
4952                                   for (iSequence = 0; iSequence < numberRows_ + numberColumns_;
4953                                             iSequence++) {
4954                                        double djValue = dj_[iSequence];
4955                                        double originalLo = originalLower(iSequence);
4956                                        double originalUp = originalUpper(iSequence);
4957                                        if (fabs(djValue) > fabs(largest)) {
4958                                             if (getStatus(iSequence) != basic) {
4959                                                  if (djValue > 0 && originalLo < -1.0e20) {
4960                                                       if (djValue > fabs(largest)) {
4961                                                            largest = djValue;
4962                                                            iChosen = iSequence;
4963                                                       }
4964                                                  } else if (djValue < 0 && originalUp > 1.0e20) {
4965                                                       if (-djValue > fabs(largest)) {
4966                                                            largest = djValue;
4967                                                            iChosen = iSequence;
4968                                                       }
4969                                                  }
4970                                             }
4971                                        }
4972                                   }
4973                              }
4974                              if (iChosen >= 0) {
4975                                   int iSave = sequenceIn_;
4976                                   sequenceIn_ = iChosen;
4977                                   unpack(rowArray_[1]);
4978                                   sequenceIn_ = iSave;
4979                                   // if dual infeasibilities then must be free vector so add in dual
4980                                   if (numberDualInfeasibilities_) {
4981                                        if (fabs(changeCost) > 1.0e-5)
4982                                          COIN_DETAIL_PRINT(printf("Odd free/unbounded combo\n"));
4983                                        changeCost += cost_[iChosen];
4984                                   }
4985                                   problemStatus_ = checkUnbounded(rowArray_[1], rowArray_[0],
4986                                                                   changeCost);
4987                                   rowArray_[1]->clear();
4988                              } else {
4989                                   problemStatus_ = -3;
4990                              }
4991                              if (problemStatus_ == 2 && perturbation_ == 101) {
4992                                   perturbation_ = 102; // stop any perturbations
4993                                   cleanDuals = 1;
4994                                   createRim4(false);
4995                                   progress_.modifyObjective(-COIN_DBL_MAX);
4996                                   problemStatus_ = -1;
4997                              }
4998                              if (problemStatus_ == 2) {
4999                                   // it is unbounded - restore solution
5000                                   // but first add in changes to non-basic
5001                                   int iColumn;
5002                                   double * original = columnArray_[0]->denseVector();
5003                                   for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
5004                                        if(getColumnStatus(iColumn) != basic)
5005                                             ray_[iColumn] +=
5006                                                  saveColumnSolution[iColumn] - original[iColumn];
5007                                        columnActivityWork_[iColumn] = original[iColumn];
5008                                   }
5009                                   CoinMemcpyN(saveRowSolution, numberRows_,
5010                                               rowActivityWork_);
5011                              }
5012                         } else {
5013                              doOriginalTolerance = 2;
5014                              rowArray_[0]->clear();
5015                         }
5016                    }
5017                    delete [] saveColumnSolution;
5018                    delete [] saveRowSolution;
5019               }
5020               if (problemStatus_ == -4 || problemStatus_ == -5) {
5021                    // may be infeasible - or may be bounds are wrong
5022                    numberChangedBounds = changeBounds(0, NULL, changeCost);
5023                    needCleanFake = true;
5024                    /* Should this be here as makes no difference to being feasible.
5025                       But seems to make a difference to run times. */
5026                    if (perturbation_ == 101 && 0) {
5027                         perturbation_ = 102; // stop any perturbations
5028                         cleanDuals = 1;
5029                         numberChangedBounds = 1;
5030                         // make sure fake bounds are back
5031                         changeBounds(1, NULL, changeCost);
5032                         needCleanFake = true;
5033                         createRim4(false);
5034                         progress_.modifyObjective(-COIN_DBL_MAX);
5035                    }
5036                    if ((numberChangedBounds <= 0 || dualBound_ > 1.0e20 ||
5037                              (largestPrimalError_ > 1.0 && dualBound_ > 1.0e17)) &&
5038                              (numberPivots < 4 || sumPrimalInfeasibilities_ > 1.0e-6)) {
5039                         problemStatus_ = 1; // infeasible
5040                         if (perturbation_ == 101) {
5041                              perturbation_ = 102; // stop any perturbations
5042                              //cleanDuals=1;
5043                              //numberChangedBounds=1;
5044                              //createRim4(false);
5045                         }
5046                         // but double check
5047                         if (!numberPrimalInfeasibilities_) {
5048                           problemStatus_=-1;
5049                           doOriginalTolerance=2;
5050                         }
5051                    } else {
5052                         problemStatus_ = -1; //iterate
5053                         cleanDuals = 1;
5054                         if (numberChangedBounds <= 0)
5055                              doOriginalTolerance = 2;
5056                         // and delete ray which has been created
5057                         delete [] ray_;
5058                         ray_ = NULL;
5059                    }
5060
5061               }
5062          } else {
5063               cleanDuals = 1;
5064          }
5065          if (problemStatus_ < 0) {
5066               if (doOriginalTolerance == 2) {
5067                    // put back original tolerance
5068                    lastCleaned = numberIterations_;
5069                    numberChanged_ = 0; // Number of variables with changed costs
5070                    handler_->message(CLP_DUAL_ORIGINAL, messages_)
5071                              << CoinMessageEol;
5072                    perturbation_ = 102; // stop any perturbations
5073#if 0
5074                    double * xcost = new double[numberRows_+numberColumns_];
5075                    double * xlower = new double[numberRows_+numberColumns_];
5076                    double * xupper = new double[numberRows_+numberColumns_];
5077                    double * xdj = new double[numberRows_+numberColumns_];
5078                    double * xsolution = new double[numberRows_+numberColumns_];
5079                    CoinMemcpyN(cost_, (numberRows_ + numberColumns_), xcost);
5080                    CoinMemcpyN(lower_, (numberRows_ + numberColumns_), xlower);
5081                    CoinMemcpyN(upper_, (numberRows_ + numberColumns_), xupper);
5082                    CoinMemcpyN(dj_, (numberRows_ + numberColumns_), xdj);
5083                    CoinMemcpyN(solution_, (numberRows_ + numberColumns_), xsolution);
5084#endif
5085                    createRim4(false);
5086                    progress_.modifyObjective(-COIN_DBL_MAX);
5087                    // make sure duals are current
5088                    computeDuals(givenDuals);
5089                    checkDualSolution();
5090#if 0
5091                    int i;
5092                    for (i = 0; i < numberRows_ + numberColumns_; i++) {
5093                         if (cost_[i] != xcost[i])
5094                              printf("** %d old cost %g new %g sol %g\n",
5095                                     i, xcost[i], cost_[i], solution_[i]);
5096                         if (lower_[i] != xlower[i])
5097                              printf("** %d old lower %g new %g sol %g\n",
5098                                     i, xlower[i], lower_[i], solution_[i]);
5099                         if (upper_[i] != xupper[i])
5100                              printf("** %d old upper %g new %g sol %g\n",
5101                                     i, xupper[i], upper_[i], solution_[i]);
5102                         if (dj_[i] != xdj[i])
5103                              printf("** %d old dj %g new %g sol %g\n",
5104                                     i, xdj[i], dj_[i], solution_[i]);
5105                         if (solution_[i] != xsolution[i])
5106                              printf("** %d old solution %g new %g sol %g\n",
5107                                     i, xsolution[i], solution_[i], solution_[i]);
5108                    }
5109                    //delete [] xcost;
5110                    //delete [] xupper;
5111                    //delete [] xlower;
5112                    //delete [] xdj;
5113                    //delete [] xsolution;
5114#endif
5115                    // put back bounds as they were if was optimal
5116                    if (doOriginalTolerance == 2 && cleanDuals != 2) {
5117                         changeMade_++; // say something changed
5118                         /* We may have already changed some bounds in this function
5119                            so save numberFake_ and add in.
5120
5121                            Worst that can happen is that we waste a bit of time  - but it must be finite.
5122                         */
5123                         //int saveNumberFake = numberFake_;
5124                         //resetFakeBounds(-1);
5125                         changeBounds(3, NULL, changeCost);
5126                         needCleanFake = true;
5127                         //numberFake_ += saveNumberFake;
5128                         //resetFakeBounds(-1);
5129                         cleanDuals = 2;
5130                         //cleanDuals=1;
5131                    }
5132#if 0
5133                    //int i;
5134                    for (i = 0; i < numberRows_ + numberColumns_; i++) {
5135                         if (cost_[i] != xcost[i])
5136                              printf("** %d old cost %g new %g sol %g\n",
5137                                     i, xcost[i], cost_[i], solution_[i]);
5138                         if (lower_[i] != xlower[i])
5139                              printf("** %d old lower %g new %g sol %g\n",
5140                                     i, xlower[i], lower_[i], solution_[i]);
5141                         if (upper_[i] != xupper[i])
5142                              printf("** %d old upper %g new %g sol %g\n",
5143                                     i, xupper[i], upper_[i], solution_[i]);
5144                         if (dj_[i] != xdj[i])
5145                              printf("** %d old dj %g new %g sol %g\n",
5146                                     i, xdj[i], dj_[i], solution_[i]);
5147                         if (solution_[i] != xsolution[i])
5148                              printf("** %d old solution %g new %g sol %g\n",
5149                                     i, xsolution[i], solution_[i], solution_[i]);
5150                    }
5151                    delete [] xcost;
5152                    delete [] xupper;
5153                    delete [] xlower;
5154                    delete [] xdj;
5155                    delete [] xsolution;
5156#endif
5157               }
5158               if (cleanDuals == 1 || (cleanDuals == 2 && !numberDualInfeasibilities_)) {
5159                    // make sure dual feasible
5160                    // look at all rows and columns
5161                    rowArray_[0]->clear();
5162                    columnArray_[0]->clear();
5163                    double objectiveChange = 0.0;
5164                    double savePrimalInfeasibilities = sumPrimalInfeasibilities_;
5165                    if (!numberIterations_) {
5166                      int nTotal = numberRows_ + numberColumns_;
5167                      if (arraysNotCreated) {
5168                        // create save arrays
5169                        delete [] saveStatus_;
5170                        delete [] savedSolution_;
5171                        saveStatus_ = new unsigned char [nTotal];
5172                        savedSolution_ = new double [nTotal];
5173                        arraysNotCreated = false;
5174                      }
5175                      // save arrays
5176                      CoinMemcpyN(status_, nTotal, saveStatus_);
5177                      CoinMemcpyN(rowActivityWork_,
5178                                  numberRows_, savedSolution_ + numberColumns_);
5179                      CoinMemcpyN(columnActivityWork_, numberColumns_, savedSolution_);
5180                    }
5181#if 0
5182                    double * xcost = new double[numberRows_+numberColumns_];
5183                    double * xlower = new double[numberRows_+numberColumns_];
5184                    double * xupper = new double[numberRows_+numberColumns_];
5185                    double * xdj = new double[numberRows_+numberColumns_];
5186                    double * xsolution = new double[numberRows_+numberColumns_];
5187                    CoinMemcpyN(cost_, (numberRows_ + numberColumns_), xcost);
5188                    CoinMemcpyN(lower_, (numberRows_ + numberColumns_), xlower);
5189                    CoinMemcpyN(upper_, (numberRows_ + numberColumns_), xupper);
5190                    CoinMemcpyN(dj_, (numberRows_ + numberColumns_), xdj);
5191                    CoinMemcpyN(solution_, (numberRows_ + numberColumns_), xsolution);
5192#endif
5193                    if (givenDuals)
5194                         dualTolerance_ = 1.0e50;
5195#if CLP_CAN_HAVE_ZERO_OBJ>1
5196                    if ((specialOptions_&2097152)==0) {
5197#endif
5198                    updateDualsInDual(rowArray_[0], columnArray_[0], rowArray_[1],
5199                                      0.0, objectiveChange, true);
5200#if CLP_CAN_HAVE_ZERO_OBJ>1
5201                    } else {
5202                      rowArray_[0]->clear();
5203                      rowArray_[1]->clear();
5204                      columnArray_[0]->clear();
5205                    }
5206#endif
5207                    dualTolerance_ = saveTolerance;
5208#if 0
5209                    int i;
5210                    for (i = 0; i < numberRows_ + numberColumns_; i++) {
5211                         if (cost_[i] != xcost[i])
5212                              printf("** %d old cost %g new %g sol %g\n",
5213                                     i, xcost[i], cost_[i], solution_[i]);
5214                         if (lower_[i] != xlower[i])
5215                              printf("** %d old lower %g new %g sol %g\n",
5216                                     i, xlower[i], lower_[i], solution_[i]);
5217                         if (upper_[i] != xupper[i])
5218                              printf("** %d old upper %g new %g sol %g\n",
5219                                     i, xupper[i], upper_[i], solution_[i]);
5220                         if (dj_[i] != xdj[i])
5221                              printf("** %d old dj %g new %g sol %g\n",
5222                                     i, xdj[i], dj_[i], solution_[i]);
5223                         if (solution_[i] != xsolution[i])
5224                              printf("** %d old solution %g new %g sol %g\n",
5225                                     i, xsolution[i], solution_[i], solution_[i]);
5226                    }
5227                    delete [] xcost;
5228                    delete [] xupper;
5229                    delete [] xlower;
5230                    delete [] xdj;
5231                    delete [] xsolution;
5232#endif
5233                    // for now - recompute all
5234                    gutsOfSolution(NULL, NULL);
5235                    if (givenDuals)
5236                         dualTolerance_ = 1.0e50;
5237#if CLP_CAN_HAVE_ZERO_OBJ>1
5238                    if ((specialOptions_&2097152)==0) {
5239#endif
5240                    updateDualsInDual(rowArray_[0], columnArray_[0], rowArray_[1],
5241                                      0.0, objectiveChange, true);
5242#if CLP_CAN_HAVE_ZERO_OBJ>1
5243                    } else {
5244                      rowArray_[0]->clear();
5245                      rowArray_[1]->clear();
5246                      columnArray_[0]->clear();