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

Last change on this file since 1924 was 1924, checked in by forrest, 7 years ago

changes for miqp, parameters and presolve

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