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

Last change on this file since 1977 was 1977, checked in by forrest, 6 years ago

fix for tiny qp

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