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

Last change on this file since 1785 was 1785, checked in by forrest, 8 years ago

patches to try and make parametrics faster

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