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

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

make faster on difficult problems

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