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

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

mainly dantzig-wolfe - also Clp strong branching (rather than OsiClp?)

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 330.4 KB
Line 
1/* $Id: ClpSimplexDual.cpp 2025 2014-03-19 12:49:55Z 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_) && (moreSpecialOptions_&1048576)==0
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_) && (moreSpecialOptions_&1048576)==0
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                                   *handler_ << "maybe forcing re-factorization - sum " << sumBadPivots << " " << factorization_->pivots() << " pivots" << CoinMessageEol;
3681                              if(factorization_->pivots() > 3) {
3682                                   badSumPivots = true;
3683                                   break;
3684                              }
3685                         }
3686#endif
3687                         swapped[1-iFlip] = numberPossiblySwapped;
3688                         interesting[1-iFlip] = numberRemaining;
3689                         // If we stop now this will be increase in objective (I think)
3690                         double increase = (fabs(dualOut_) - totalThru) * theta_;
3691                         increase += increaseInObjective;
3692                         if (theta_ < 0.0)
3693                              thruThis += fabs(dualOut_); // force using this one
3694                         if (increaseInObjective < 0.0 && increase < 0.0 && lastSequence >= 0) {
3695                              // back
3696                              // We may need to be more careful - we could do by
3697                              // switch so we always do fine grained?
3698                              bestPivot = 0.0;
3699                         } else {
3700                              // add in
3701                              totalThru += thruThis;
3702                              increaseInObjective += increaseInThis;
3703                         }
3704                         if (bestPivot < 0.1 * bestEverPivot &&
3705                                   bestEverPivot > 1.0e-6 &&
3706                                   (bestPivot < 1.0e-3 || totalThru * 2.0 > fabs(dualOut_))) {
3707                              // back to previous one
3708                              sequenceIn_ = lastSequence;
3709                              // swap regions
3710                              iFlip = 1 - iFlip;
3711                              break;
3712                         } else if (sequenceIn_ == -1 && upperTheta > largeValue_) {
3713                              if (lastPivot > acceptablePivot) {
3714                                   // back to previous one
3715                                   sequenceIn_ = lastSequence;
3716                                   // swap regions
3717                                   iFlip = 1 - iFlip;
3718                              } else {
3719                                   // can only get here if all pivots too small
3720                              }
3721                              break;
3722                         } else if (totalThru >= fabs(dualOut_)) {
3723                              modifyCosts = true; // fine grain - we can modify costs
3724                              break; // no point trying another loop
3725                         } else {
3726                              lastSequence = sequenceIn_;
3727                              if (bestPivot > bestEverPivot)
3728                                   bestEverPivot = bestPivot;
3729                              iFlip = 1 - iFlip;
3730                              modifyCosts = true; // fine grain - we can modify costs
3731                         }
3732                    }
3733                    if (iTry == MAXTRY)
3734                         iFlip = 1 - iFlip; // flip back
3735                    break;
3736               } else {
3737                    // skip this lot
3738                    if (bestPivot > 1.0e-3 || bestPivot > bestEverPivot) {
3739                         bestEverPivot = bestPivot;
3740                         lastSequence = bestSequence;
3741                    } else {
3742                         // keep old swapped
3743                         CoinMemcpyN(array[iFlip] + swapped[iFlip],
3744                                     numberColumns_ - swapped[iFlip], array[1-iFlip] + swapped[iFlip]);
3745                         CoinMemcpyN(indices[iFlip] + swapped[iFlip],
3746                                     numberColumns_ - swapped[iFlip], indices[1-iFlip] + swapped[iFlip]);
3747                         marker[1-iFlip][1] = CoinMin(marker[1-iFlip][1], swapped[iFlip]);
3748                         swapped[1-iFlip] = swapped[iFlip];
3749                    }
3750                    increaseInObjective += increaseInThis;
3751                    iFlip = 1 - iFlip; // swap regions
3752                    tentativeTheta = 2.0 * upperTheta;
3753                    totalThru += thruThis;
3754               }
3755          }
3756
3757          // can get here without sequenceIn_ set but with lastSequence
3758          if (sequenceIn_ < 0 && lastSequence >= 0) {
3759               // back to previous one
3760               sequenceIn_ = lastSequence;
3761               // swap regions
3762               iFlip = 1 - iFlip;
3763          }
3764
3765#define MINIMUMTHETA 1.0e-18
3766          // Movement should be minimum for anti-degeneracy - unless
3767          // fixed variable out
3768          double minimumTheta;
3769          if (upperOut_ > lowerOut_)
3770               minimumTheta = MINIMUMTHETA;
3771          else
3772               minimumTheta = 0.0;
3773          if (sequenceIn_ >= 0) {
3774               // at this stage sequenceIn_ is just pointer into index array
3775               // flip just so we can use iFlip
3776               iFlip = 1 - iFlip;
3777               spare = array[iFlip];
3778               index = indices[iFlip];
3779               double oldValue;
3780               alpha_ = spare[sequenceIn_];
3781               sequenceIn_ = indices[iFlip][sequenceIn_];
3782               oldValue = dj_[sequenceIn_];
3783               theta_ = CoinMax(oldValue / alpha_, 0.0);
3784               if (theta_ < minimumTheta && fabs(alpha_) < 1.0e5 && 1) {
3785                    // can't pivot to zero
3786#if 0
3787                    if (oldValue - minimumTheta*alpha_ >= -dualTolerance_) {
3788                         theta_ = minimumTheta;
3789                    } else if (oldValue - minimumTheta*alpha_ >= -newTolerance) {
3790                         theta_ = minimumTheta;
3791                         thisIncrease = true;
3792                    } else {
3793                         theta_ = CoinMax((oldValue + newTolerance) / alpha_, 0.0);
3794                         thisIncrease = true;
3795                    }
3796#else
3797                    theta_ = minimumTheta;
3798#endif
3799               }
3800               // may need to adjust costs so all dual feasible AND pivoted is exactly 0
3801               //int costOffset = numberRows_+numberColumns_;
3802               if (modifyCosts && !badSumPivots) {
3803                    int i;
3804                    for (i = numberColumns_ - 1; i >= swapped[iFlip]; i--) {
3805                         int iSequence = index[i];
3806                         double alpha = spare[i];
3807                         double value = dj_[iSequence] - theta_ * alpha;
3808
3809                         // can't be free here
3810
3811                         if (alpha < 0.0) {
3812                              //at upper bound
3813                              if (value > dualTolerance_) {
3814                                   //thisIncrease = true;
3815#if CLP_CAN_HAVE_ZERO_OBJ<2
3816#define MODIFYCOST 2
3817#endif
3818#if MODIFYCOST
3819                                   // modify cost to hit new tolerance
3820                                   double modification = alpha * theta_ - dj_[iSequence]
3821                                                         + newTolerance;
3822                                   if ((specialOptions_&(2048 + 4096 + 16384)) != 0) {
3823                                        if ((specialOptions_ & 16384) != 0) {
3824                                             if (fabs(modification) < 1.0e-8)
3825                                                  modification = 0.0;
3826                                        } else if ((specialOptions_ & 2048) != 0) {
3827                                             if (fabs(modification) < 1.0e-10)
3828                                                  modification = 0.0;
3829                                        } else {
3830                                             if (fabs(modification) < 1.0e-12)
3831                                                  modification = 0.0;
3832                                        }
3833                                   }
3834                                   dj_[iSequence] += modification;
3835                                   cost_[iSequence] +=  modification;
3836                                   if (modification)
3837                                        numberChanged_ ++; // Say changed costs
3838                                   //cost_[iSequence+costOffset] += modification; // save change
3839#endif
3840                              }
3841                         } else {
3842                              // at lower bound
3843                              if (-value > dualTolerance_) {
3844                                   //thisIncrease = true;
3845#if MODIFYCOST
3846                                   // modify cost to hit new tolerance
3847                                   double modification = alpha * theta_ - dj_[iSequence]
3848                                                         - newTolerance;
3849                                   //modification = CoinMax(modification,-dualTolerance_);
3850                                   //assert (fabs(modification)<1.0e-7);
3851                                   if ((specialOptions_&(2048 + 4096)) != 0) {
3852                                        if ((specialOptions_ & 2048) != 0) {
3853                                             if (fabs(modification) < 1.0e-10)
3854                                                  modification = 0.0;
3855                                        } else {
3856                                             if (fabs(modification) < 1.0e-12)
3857                                                  modification = 0.0;
3858                                        }
3859                                   }
3860                                   dj_[iSequence] += modification;
3861                                   cost_[iSequence] +=  modification;
3862                                   if (modification)
3863                                        numberChanged_ ++; // Say changed costs
3864                                   //cost_[iSequence+costOffset] += modification; // save change
3865#endif
3866                              }
3867                         }
3868                    }
3869               }
3870          }
3871     }
3872
3873#ifdef MORE_CAREFUL
3874     // If we have done pivots and things look bad set alpha_ 0.0 to force factorization
3875     if ((badSumPivots ||
3876               fabs(theta_ * badFree) > 10.0 * dualTolerance_) && factorization_->pivots()) {
3877          if (handler_->logLevel() > 1)
3878               *handler_ << "forcing re-factorization" << CoinMessageEol;
3879          sequenceIn_ = -1;
3880     }
3881#endif
3882     if (sequenceIn_ >= 0) {
3883          lowerIn_ = lower_[sequenceIn_];
3884          upperIn_ = upper_[sequenceIn_];
3885          valueIn_ = solution_[sequenceIn_];
3886          dualIn_ = dj_[sequenceIn_];
3887
3888          if (numberTimesOptimal_) {
3889               // can we adjust cost back closer to original
3890               //*** add coding
3891          }
3892#if MODIFYCOST>1
3893          // modify cost to hit zero exactly
3894          // so (dualIn_+modification)==theta_*alpha_
3895          double modification = theta_ * alpha_ - dualIn_;
3896          // But should not move objective too much ??
3897#define DONT_MOVE_OBJECTIVE
3898#ifdef DONT_MOVE_OBJECTIVE
3899          double moveObjective = fabs(modification * solution_[sequenceIn_]);
3900          double smallMove = CoinMax(fabs(objectiveValue_), 1.0e-3);
3901          if (moveObjective > smallMove) {
3902               if (handler_->logLevel() > 1)
3903                    printf("would move objective by %g - original mod %g sol value %g\n", moveObjective,
3904                           modification, solution_[sequenceIn_]);
3905               modification *= smallMove / moveObjective;
3906          }
3907#endif
3908          if (badSumPivots)
3909               modification = 0.0;
3910          if ((specialOptions_&(2048 + 4096)) != 0) {
3911               if ((specialOptions_ & 16384) != 0) {
3912                    // in fast dual
3913                    if (fabs(modification) < 1.0e-7)
3914                         modification = 0.0;
3915               } else if ((specialOptions_ & 2048) != 0) {
3916                    if (fabs(modification) < 1.0e-10)
3917                         modification = 0.0;
3918               } else {
3919                    if (fabs(modification) < 1.0e-12)
3920                         modification = 0.0;
3921               }
3922          }
3923          dualIn_ += modification;
3924          dj_[sequenceIn_] = dualIn_;
3925          cost_[sequenceIn_] += modification;
3926          if (modification)
3927               numberChanged_ ++; // Say changed costs
3928          //int costOffset = numberRows_+numberColumns_;
3929          //cost_[sequenceIn_+costOffset] += modification; // save change
3930          //assert (fabs(modification)<1.0e-6);
3931#ifdef CLP_DEBUG
3932          if ((handler_->logLevel() & 32) && fabs(modification) > 1.0e-15)
3933               printf("exact %d new cost %g, change %g\n", sequenceIn_,
3934                      cost_[sequenceIn_], modification);
3935#endif
3936#endif
3937
3938          if (alpha_ < 0.0) {
3939               // as if from upper bound
3940               directionIn_ = -1;
3941               upperIn_ = valueIn_;
3942          } else {
3943               // as if from lower bound
3944               directionIn_ = 1;
3945               lowerIn_ = valueIn_;
3946          }
3947     } else {
3948          // no pivot
3949          bestPossible = 0.0;
3950          alpha_ = 0.0;
3951     }
3952     //if (thisIncrease)
3953     //dualTolerance_+= 1.0e-6*dblParam_[ClpDualTolerance];
3954
3955     // clear arrays
3956
3957     for (i = 0; i < 2; i++) {
3958          CoinZeroN(array[i], marker[i][0]);
3959          CoinZeroN(array[i] + marker[i][1], numberColumns_ - marker[i][1]);
3960     }
3961     return bestPossible;
3962}
3963#ifdef CLP_ALL_ONE_FILE
3964#undef MAXTRY
3965#endif
3966/* Checks if tentative optimal actually means unbounded
3967   Returns -3 if not, 2 if is unbounded */
3968int
3969ClpSimplexDual::checkUnbounded(CoinIndexedVector * ray,
3970                               CoinIndexedVector * spare,
3971                               double changeCost)
3972{
3973     int status = 2; // say unbounded
3974     factorization_->updateColumn(spare, ray);
3975     // get reduced cost
3976     int i;
3977     int number = ray->getNumElements();
3978     int * index = ray->getIndices();
3979     double * array = ray->denseVector();
3980     for (i = 0; i < number; i++) {
3981          int iRow = index[i];
3982          int iPivot = pivotVariable_[iRow];
3983          changeCost -= cost(iPivot) * array[iRow];
3984     }
3985     double way;
3986     if (changeCost > 0.0) {
3987          //try going down
3988          way = 1.0;
3989     } else if (changeCost < 0.0) {
3990          //try going up
3991          way = -1.0;
3992     } else {
3993#ifdef CLP_DEBUG
3994          printf("can't decide on up or down\n");
3995#endif
3996          way = 0.0;
3997          status = -3;
3998     }
3999     double movement = 1.0e10 * way; // some largish number
4000     double zeroTolerance = 1.0e-14 * dualBound_;
4001     for (i = 0; i < number; i++) {
4002          int iRow = index[i];
4003          int iPivot = pivotVariable_[iRow];
4004          double arrayValue = array[iRow];
4005          if (fabs(arrayValue) < zeroTolerance)
4006               arrayValue = 0.0;
4007          double newValue = solution(iPivot) + movement * arrayValue;
4008          if (newValue > upper(iPivot) + primalTolerance_ ||
4009                    newValue < lower(iPivot) - primalTolerance_)
4010               status = -3; // not unbounded
4011     }
4012     if (status == 2) {
4013          // create ray
4014          delete [] ray_;
4015          ray_ = new double [numberColumns_];
4016          CoinZeroN(ray_, numberColumns_);
4017          for (i = 0; i < number; i++) {
4018               int iRow = index[i];
4019               int iPivot = pivotVariable_[iRow];
4020               double arrayValue = array[iRow];
4021               if (iPivot < numberColumns_ && fabs(arrayValue) >= zeroTolerance)
4022                    ray_[iPivot] = way * array[iRow];
4023          }
4024     }
4025     ray->clear();
4026     return status;
4027}
4028//static int count_alpha=0;
4029/* Checks if finished.  Updates status */
4030void
4031ClpSimplexDual::statusOfProblemInDual(int & lastCleaned, int type,
4032                                      double * givenDuals, ClpDataSave & saveData,
4033                                      int ifValuesPass)
4034{
4035#ifdef CLP_INVESTIGATE_SERIAL
4036     if (z_thinks > 0 && z_thinks < 2)
4037          z_thinks += 2;
4038#endif
4039     bool arraysNotCreated = (type==0);
4040     // If lots of iterations then adjust costs if large ones
4041     if (numberIterations_ > 4 * (numberRows_ + numberColumns_) && objectiveScale_ == 1.0) {
4042          double largest = 0.0;
4043          for (int i = 0; i < numberRows_; i++) {
4044               int iColumn = pivotVariable_[i];
4045               largest = CoinMax(largest, fabs(cost_[iColumn]));
4046          }
4047          if (largest > 1.0e6) {
4048               objectiveScale_ = 1.0e6 / largest;
4049               for (int i = 0; i < numberRows_ + numberColumns_; i++)
4050                    cost_[i] *= objectiveScale_;
4051          }
4052     }
4053     int numberPivots = factorization_->pivots();
4054     double realDualInfeasibilities = 0.0;
4055     if (type == 2) {
4056          if (alphaAccuracy_ != -1.0)
4057               alphaAccuracy_ = -2.0;
4058          // trouble - restore solution
4059          CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4060          CoinMemcpyN(savedSolution_ + numberColumns_ ,
4061                      numberRows_, rowActivityWork_);
4062          CoinMemcpyN(savedSolution_ ,
4063                      numberColumns_, columnActivityWork_);
4064          // restore extra stuff
4065          int dummy;
4066          matrix_->generalExpanded(this, 6, dummy);
4067          forceFactorization_ = 1; // a bit drastic but ..
4068          changeMade_++; // say something changed
4069          // get correct bounds on all variables
4070          resetFakeBounds(0);
4071     }
4072     int tentativeStatus = problemStatus_;
4073     double changeCost;
4074     bool unflagVariables = true;
4075     bool weightsSaved = false;
4076     bool weightsSaved2 = numberIterations_ && !numberPrimalInfeasibilities_;
4077     int dontFactorizePivots = dontFactorizePivots_;
4078     if (type == 3) {
4079          type = 1;
4080          dontFactorizePivots = 1;
4081     }
4082     if (alphaAccuracy_ < 0.0 || !numberPivots || alphaAccuracy_ > 1.0e4 || numberPivots > 20) {
4083          if (problemStatus_ > -3 || numberPivots > dontFactorizePivots) {
4084               // factorize
4085               // later on we will need to recover from singularities
4086               // also we could skip if first time
4087               // save dual weights
4088               dualRowPivot_->saveWeights(this, 1);
4089               weightsSaved = true;
4090               if (type) {
4091                    // is factorization okay?
4092                    if (internalFactorize(1)) {
4093                         // no - restore previous basis
4094                         unflagVariables = false;
4095                         assert (type == 1);
4096                         changeMade_++; // say something changed
4097                         // Keep any flagged variables
4098                         int i;
4099                         for (i = 0; i < numberRows_ + numberColumns_; i++) {
4100                              if (flagged(i))
4101                                   saveStatus_[i] |= 64; //say flagged
4102                         }
4103                         CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4104                         CoinMemcpyN(savedSolution_ + numberColumns_ ,
4105                                     numberRows_, rowActivityWork_);
4106                         CoinMemcpyN(savedSolution_ ,
4107                                     numberColumns_, columnActivityWork_);
4108                         // restore extra stuff
4109                         int dummy;
4110                         matrix_->generalExpanded(this, 6, dummy);
4111                         // get correct bounds on all variables
4112                         resetFakeBounds(1);
4113                         // need to reject something
4114                         char x = isColumn(sequenceOut_) ? 'C' : 'R';
4115                         handler_->message(CLP_SIMPLEX_FLAG, messages_)
4116                                   << x << sequenceWithin(sequenceOut_)
4117                                   << CoinMessageEol;
4118#ifdef COIN_DEVELOP
4119                         printf("flag d\n");
4120#endif
4121                         setFlagged(sequenceOut_);
4122                         progress_.clearBadTimes();
4123
4124                         // Go to safe
4125                         factorization_->pivotTolerance(0.99);
4126                         forceFactorization_ = 1; // a bit drastic but ..
4127                         type = 2;
4128                         //assert (internalFactorize(1)==0);
4129                         if (internalFactorize(1)) {
4130                              CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4131                              CoinMemcpyN(savedSolution_ + numberColumns_ ,
4132                                          numberRows_, rowActivityWork_);
4133                              CoinMemcpyN(savedSolution_ ,
4134                                          numberColumns_, columnActivityWork_);
4135                              // restore extra stuff
4136                              int dummy;
4137                              matrix_->generalExpanded(this, 6, dummy);
4138                              // debug
4139                              int returnCode = internalFactorize(1);
4140                              while (returnCode) {
4141                                   // ouch
4142                                   // switch off dense
4143                                   int saveDense = factorization_->denseThreshold();
4144                                   factorization_->setDenseThreshold(0);
4145                                   // Go to safe
4146                                   factorization_->pivotTolerance(0.99);
4147                                   // make sure will do safe factorization
4148                                   pivotVariable_[0] = -1;
4149                                   returnCode = internalFactorize(2);
4150                                   factorization_->setDenseThreshold(saveDense);
4151                              }
4152                              // get correct bounds on all variables
4153                              resetFakeBounds(1);
4154                         }
4155                    }
4156               }
4157               if (problemStatus_ != -4 || numberPivots > 10)
4158                    problemStatus_ = -3;
4159          }
4160     } else {
4161          //printf("testing with accuracy of %g and status of %d\n",alphaAccuracy_,problemStatus_);
4162          //count_alpha++;
4163          //if ((count_alpha%5000)==0)
4164          //printf("count alpha %d\n",count_alpha);
4165     }
4166     // at this stage status is -3 or -4 if looks infeasible
4167     // get primal and dual solutions
4168#if 0
4169     {
4170          int numberTotal = numberRows_ + numberColumns_;
4171          double * saveSol = CoinCopyOfArray(solution_, numberTotal);
4172          double * saveDj = CoinCopyOfArray(dj_, numberTotal);
4173          double tolerance = type ? 1.0e-4 : 1.0e-8;
4174          // always if values pass
4175          double saveObj = objectiveValue_;
4176          double sumPrimal = sumPrimalInfeasibilities_;
4177          int numberPrimal = numberPrimalInfeasibilities_;
4178          double sumDual = sumDualInfeasibilities_;
4179          int numberDual = numberDualInfeasibilities_;
4180          gutsOfSolution(givenDuals, NULL);
4181          int j;
4182          double largestPrimal = tolerance;
4183          int iPrimal = -1;
4184          for (j = 0; j < numberTotal; j++) {
4185               double difference = solution_[j] - saveSol[j];
4186               if (fabs(difference) > largestPrimal) {
4187                    iPrimal = j;
4188                    largestPrimal = fabs(difference);
4189               }
4190          }
4191          double largestDual = tolerance;
4192          int iDual = -1;
4193          for (j = 0; j < numberTotal; j++) {
4194               double difference = dj_[j] - saveDj[j];
4195               if (fabs(difference) > largestDual && upper_[j] > lower_[j]) {
4196                    iDual = j;
4197                    largestDual = fabs(difference);
4198               }
4199          }
4200          if (!type) {
4201               if (fabs(saveObj - objectiveValue_) > 1.0e-5 ||
4202                         numberPrimal != numberPrimalInfeasibilities_ || numberPrimal != 1 ||
4203                         fabs(sumPrimal - sumPrimalInfeasibilities_) > 1.0e-5 || iPrimal >= 0 ||
4204                         numberDual != numberDualInfeasibilities_ || numberDual != 0 ||
4205                         fabs(sumDual - sumDualInfeasibilities_) > 1.0e-5 || iDual >= 0)
4206                    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",
4207                           type, numberIterations_, numberPivots,
4208                           numberPrimal, numberPrimalInfeasibilities_, sumPrimal, sumPrimalInfeasibilities_,
4209                           largestPrimal, iPrimal,
4210                           numberDual, numberDualInfeasibilities_, sumDual, sumDualInfeasibilities_,
4211                           largestDual, iDual,
4212                           saveObj, objectiveValue_);
4213          } else {
4214               if (fabs(saveObj - objectiveValue_) > 1.0e-5 ||
4215                         numberPrimalInfeasibilities_ || iPrimal >= 0 ||
4216                         numberDualInfeasibilities_ || iDual >= 0)
4217                    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",
4218                           type, numberIterations_, numberPivots,
4219                           numberPrimal, numberPrimalInfeasibilities_, sumPrimal, sumPrimalInfeasibilities_,
4220                           largestPrimal, iPrimal,
4221                           numberDual, numberDualInfeasibilities_, sumDual, sumDualInfeasibilities_,
4222                           largestDual, iDual,
4223                           saveObj, objectiveValue_);
4224          }
4225          delete [] saveSol;
4226          delete [] saveDj;
4227     }
4228#else
4229     if (type || ifValuesPass)
4230          gutsOfSolution(givenDuals, NULL);
4231#endif
4232     // If bad accuracy treat as singular
4233     if ((largestPrimalError_ > 1.0e15 || largestDualError_ > 1.0e15) && numberIterations_) {
4234          // restore previous basis
4235          unflagVariables = false;
4236          changeMade_++; // say something changed
4237          // Keep any flagged variables
4238          int i;
4239          for (i = 0; i < numberRows_ + numberColumns_; i++) {
4240               if (flagged(i))
4241                    saveStatus_[i] |= 64; //say flagged
4242          }
4243          CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4244          CoinMemcpyN(savedSolution_ + numberColumns_ ,
4245                      numberRows_, rowActivityWork_);
4246          CoinMemcpyN(savedSolution_ ,
4247                      numberColumns_, columnActivityWork_);
4248          // restore extra stuff
4249          int dummy;
4250          matrix_->generalExpanded(this, 6, dummy);
4251          // get correct bounds on all variables
4252          resetFakeBounds(1);
4253          // need to reject something
4254          char x = isColumn(sequenceOut_) ? 'C' : 'R';
4255          handler_->message(CLP_SIMPLEX_FLAG, messages_)
4256                    << x << sequenceWithin(sequenceOut_)
4257                    << CoinMessageEol;
4258#ifdef COIN_DEVELOP
4259          printf("flag e\n");
4260#endif
4261          setFlagged(sequenceOut_);
4262          progress_.clearBadTimes();
4263
4264          // Go to safer
4265          double newTolerance = CoinMin(1.1 * factorization_->pivotTolerance(), 0.99);
4266          factorization_->pivotTolerance(newTolerance);
4267          forceFactorization_ = 1; // a bit drastic but ..
4268          if (alphaAccuracy_ != -1.0)
4269               alphaAccuracy_ = -2.0;
4270          type = 2;
4271          //assert (internalFactorize(1)==0);
4272          if (internalFactorize(1)) {
4273               CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4274               CoinMemcpyN(savedSolution_ + numberColumns_ ,
4275                           numberRows_, rowActivityWork_);
4276               CoinMemcpyN(savedSolution_ ,
4277                           numberColumns_, columnActivityWork_);
4278               // restore extra stuff
4279               int dummy;
4280               matrix_->generalExpanded(this, 6, dummy);
4281               // debug
4282               int returnCode = internalFactorize(1);
4283               while (returnCode) {
4284                    // ouch
4285                    // switch off dense
4286                    int saveDense = factorization_->denseThreshold();
4287                    factorization_->setDenseThreshold(0);
4288                    // Go to safe
4289                    factorization_->pivotTolerance(0.99);
4290                    // make sure will do safe factorization
4291                    pivotVariable_[0] = -1;
4292                    returnCode = internalFactorize(2);
4293                    factorization_->setDenseThreshold(saveDense);
4294               }
4295               // get correct bounds on all variables
4296               resetFakeBounds(1);
4297          }
4298          // get primal and dual solutions
4299          gutsOfSolution(givenDuals, NULL);
4300     } else if (goodAccuracy()) {
4301          // Can reduce tolerance
4302          double newTolerance = CoinMax(0.99 * factorization_->pivotTolerance(), saveData.pivotTolerance_);
4303          factorization_->pivotTolerance(newTolerance);
4304     }
4305     bestObjectiveValue_ = CoinMax(bestObjectiveValue_,
4306                                   objectiveValue_ - bestPossibleImprovement_);
4307     bool reallyBadProblems = false;
4308     // Double check infeasibility if no action
4309     if (progress_.lastIterationNumber(0) == numberIterations_) {
4310          if (dualRowPivot_->looksOptimal()) {
4311               numberPrimalInfeasibilities_ = 0;
4312               sumPrimalInfeasibilities_ = 0.0;
4313          }
4314#if 1
4315     } else {
4316          double thisObj = objectiveValue_ - bestPossibleImprovement_;
4317#ifdef CLP_INVESTIGATE
4318          assert (bestPossibleImprovement_ > -1000.0 && objectiveValue_ > -1.0e100);
4319          if (bestPossibleImprovement_)
4320               printf("obj %g add in %g -> %g\n", objectiveValue_, bestPossibleImprovement_,
4321                      thisObj);
4322#endif
4323          double lastObj = progress_.lastObjective(0);
4324#ifndef NDEBUG
4325#ifdef COIN_DEVELOP
4326          resetFakeBounds(-1);
4327#endif
4328#endif
4329#ifdef CLP_REPORT_PROGRESS
4330          ixxxxxx++;
4331          if (ixxxxxx >= ixxyyyy - 4 && ixxxxxx <= ixxyyyy) {
4332               char temp[20];
4333               sprintf(temp, "sol%d.out", ixxxxxx);
4334               printf("sol%d.out\n", ixxxxxx);
4335               FILE * fp = fopen(temp, "w");
4336               int nTotal = numberRows_ + numberColumns_;
4337               for (int i = 0; i < nTotal; i++)
4338                    fprintf(fp, "%d %d %g %g %g %g %g\n",
4339                            i, status_[i], lower_[i], solution_[i], upper_[i], cost_[i], dj_[i]);
4340               fclose(fp);
4341          }
4342#endif
4343          if(!ifValuesPass && firstFree_ < 0) {
4344               double testTol = 5.0e-3;
4345               if (progress_.timesFlagged() > 10) {
4346                    testTol *= pow(2.0, progress_.timesFlagged() - 8);
4347               } else if (progress_.timesFlagged() > 5) {
4348                    testTol *= 5.0;
4349               }
4350               if (lastObj > thisObj +
4351                         testTol*(fabs(thisObj) + fabs(lastObj)) + testTol) {
4352                    int maxFactor = factorization_->maximumPivots();
4353                    if ((specialOptions_ & 1048576) == 0) {
4354                         if (progress_.timesFlagged() > 10)
4355                              progress_.incrementReallyBadTimes();
4356                         if (maxFactor > 10 - 9) {
4357#ifdef COIN_DEVELOP
4358                              printf("lastobj %g thisobj %g\n", lastObj, thisObj);
4359#endif
4360                              //if (forceFactorization_<0)
4361                              //forceFactorization_= maxFactor;
4362                              //forceFactorization_ = CoinMax(1,(forceFactorization_>>1));
4363                              if ((progressFlag_ & 4) == 0 && lastObj < thisObj + 1.0e4 &&
4364                                        largestPrimalError_ < 1.0e2) {
4365                                   // Just save costs
4366                                   // save extra copy of cost_
4367                                   int nTotal = numberRows_ + numberColumns_;
4368                                   double * temp = new double [2*nTotal];
4369                                   memcpy(temp, cost_, nTotal * sizeof(double));
4370                                   memcpy(temp + nTotal, cost_, nTotal * sizeof(double));
4371                                   delete [] cost_;
4372                                   cost_ = temp;
4373                                   objectiveWork_ = cost_;
4374                                   rowObjectiveWork_ = cost_ + numberColumns_;
4375                                   progressFlag_ |= 4;
4376                              } else {
4377                                   forceFactorization_ = 1;
4378#ifdef COIN_DEVELOP
4379                                   printf("Reducing factorization frequency - bad backwards\n");
4380#endif
4381#if 1
4382                                   unflagVariables = false;
4383                                   changeMade_++; // say something changed
4384                                   int nTotal = numberRows_ + numberColumns_;
4385                                   CoinMemcpyN(saveStatus_, nTotal, status_);
4386                                   CoinMemcpyN(savedSolution_ + numberColumns_ ,
4387                                               numberRows_, rowActivityWork_);
4388                                   CoinMemcpyN(savedSolution_ ,
4389                                               numberColumns_, columnActivityWork_);
4390                                   if ((progressFlag_ & 4) == 0) {
4391                                        // save extra copy of cost_
4392                                        double * temp = new double [2*nTotal];
4393                                        memcpy(temp, cost_, nTotal * sizeof(double));
4394                                        memcpy(temp + nTotal, cost_, nTotal * sizeof(double));
4395                                        delete [] cost_;
4396                                        cost_ = temp;
4397                                        objectiveWork_ = cost_;
4398                                        rowObjectiveWork_ = cost_ + numberColumns_;
4399                                        progressFlag_ |= 4;
4400                                   } else {
4401                                        memcpy(cost_, cost_ + nTotal, nTotal * sizeof(double));
4402                                   }
4403                                   // restore extra stuff
4404                                   int dummy;
4405                                   matrix_->generalExpanded(this, 6, dummy);
4406                                   double pivotTolerance = factorization_->pivotTolerance();
4407                                   if(pivotTolerance < 0.2)
4408                                        factorization_->pivotTolerance(0.2);
4409                                   else if(progress_.timesFlagged() > 2)
4410                                        factorization_->pivotTolerance(CoinMin(pivotTolerance * 1.1, 0.99));
4411                                   if (alphaAccuracy_ != -1.0)
4412                                        alphaAccuracy_ = -2.0;
4413                                   if (internalFactorize(1)) {
4414                                        CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4415                                        CoinMemcpyN(savedSolution_ + numberColumns_ ,
4416                                                    numberRows_, rowActivityWork_);
4417                                        CoinMemcpyN(savedSolution_ ,
4418                                                    numberColumns_, columnActivityWork_);
4419                                        // restore extra stuff
4420                                        int dummy;
4421                                        matrix_->generalExpanded(this, 6, dummy);
4422                                        // debug
4423                                        int returnCode = internalFactorize(1);
4424                                        while (returnCode) {
4425                                             // ouch
4426                                             // switch off dense
4427                                             int saveDense = factorization_->denseThreshold();
4428                                             factorization_->setDenseThreshold(0);
4429                                             // Go to safe
4430                                             factorization_->pivotTolerance(0.99);
4431                                             // make sure will do safe factorization
4432                                             pivotVariable_[0] = -1;
4433                                             returnCode = internalFactorize(2);
4434                                             factorization_->setDenseThreshold(saveDense);
4435                                        }
4436                                   }
4437                                   resetFakeBounds(0);
4438                                   type = 2; // so will restore weights
4439                                   // get primal and dual solutions
4440                                   gutsOfSolution(givenDuals, NULL);
4441                                   if (numberPivots < 2) {
4442                                        // need to reject something
4443                                        char x = isColumn(sequenceOut_) ? 'C' : 'R';
4444                                        handler_->message(CLP_SIMPLEX_FLAG, messages_)
4445                                                  << x << sequenceWithin(sequenceOut_)
4446                                                  << CoinMessageEol;
4447#ifdef COIN_DEVELOP
4448                                        printf("flag d\n");
4449#endif
4450                                        setFlagged(sequenceOut_);
4451                                        progress_.clearBadTimes();
4452                                        progress_.incrementTimesFlagged();
4453                                   }
4454                                   if (numberPivots < 10)
4455                                        reallyBadProblems = true;
4456#ifdef COIN_DEVELOP
4457                                   printf("obj now %g\n", objectiveValue_);
4458#endif
4459                                   progress_.modifyObjective(objectiveValue_
4460                                                             - bestPossibleImprovement_);
4461#endif
4462                              }
4463                         }
4464                    } else {
4465                         // in fast dual give up
4466#ifdef COIN_DEVELOP
4467                         printf("In fast dual?\n");
4468#endif
4469                         problemStatus_ = 3;
4470                    }
4471               } else if (lastObj < thisObj - 1.0e-5 * CoinMax(fabs(thisObj), fabs(lastObj)) - 1.0e-3) {
4472                    numberTimesOptimal_ = 0;
4473               }
4474          }
4475#endif
4476     }
4477     // Up tolerance if looks a bit odd
4478     if (numberIterations_ > CoinMax(1000, numberRows_ >> 4) && (specialOptions_ & 64) != 0) {
4479          if (sumPrimalInfeasibilities_ && sumPrimalInfeasibilities_ < 1.0e5) {
4480               int backIteration = progress_.lastIterationNumber(CLP_PROGRESS - 1);
4481               if (backIteration > 0 && numberIterations_ - backIteration < 9 * CLP_PROGRESS) {
4482                    if (factorization_->pivotTolerance() < 0.9) {
4483                         // up tolerance
4484                         factorization_->pivotTolerance(CoinMin(factorization_->pivotTolerance() * 1.05 + 0.02, 0.91));
4485                         //printf("tol now %g\n",factorization_->pivotTolerance());
4486                         progress_.clearIterationNumbers();
4487                    }
4488               }
4489          }
4490     }
4491     // Check if looping
4492     int loop;
4493     if (!givenDuals && type != 2)
4494          loop = progress_.looping();
4495     else
4496          loop = -1;
4497     if (progress_.reallyBadTimes() > 10) {
4498          problemStatus_ = 10; // instead - try other algorithm
4499#if COIN_DEVELOP>2
4500          printf("returning at %d\n", __LINE__);
4501#endif
4502     }
4503     int situationChanged = 0;
4504     if (loop >= 0) {
4505          problemStatus_ = loop; //exit if in loop
4506          if (!problemStatus_) {
4507               // declaring victory
4508               numberPrimalInfeasibilities_ = 0;
4509               sumPrimalInfeasibilities_ = 0.0;
4510          } else {
4511               problemStatus_ = 10; // instead - try other algorithm
4512#if COIN_DEVELOP>2
4513               printf("returning at %d\n", __LINE__);
4514#endif
4515          }
4516          return;
4517     } else if (loop < -1) {
4518          // something may have changed
4519          gutsOfSolution(NULL, NULL);
4520          situationChanged = 1;
4521     }
4522     // really for free variables in
4523     if((progressFlag_ & 2) != 0) {
4524          situationChanged = 2;
4525     }
4526     progressFlag_ &= (~3); //reset progress flag
4527     if ((progressFlag_ & 4) != 0) {
4528          // save copy of cost_
4529          int nTotal = numberRows_ + numberColumns_;
4530          memcpy(cost_ + nTotal, cost_, nTotal * sizeof(double));
4531     }
4532     /*if (!numberIterations_&&sumDualInfeasibilities_)
4533       printf("OBJ %g sumPinf %g sumDinf %g\n",
4534        objectiveValue(),sumPrimalInfeasibilities_,
4535        sumDualInfeasibilities_);*/
4536     // mark as having gone optimal if looks like it
4537     if (!numberPrimalInfeasibilities_&&
4538         !numberDualInfeasibilities_)
4539       progressFlag_ |= 8;
4540     if (handler_->detail(CLP_SIMPLEX_STATUS, messages_) < 100) {
4541          handler_->message(CLP_SIMPLEX_STATUS, messages_)
4542                    << numberIterations_ << objectiveValue();
4543          handler_->printing(sumPrimalInfeasibilities_ > 0.0)
4544                    << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
4545          handler_->printing(sumDualInfeasibilities_ > 0.0)
4546                    << sumDualInfeasibilities_ << numberDualInfeasibilities_;
4547          handler_->printing(numberDualInfeasibilitiesWithoutFree_
4548                             < numberDualInfeasibilities_)
4549                    << numberDualInfeasibilitiesWithoutFree_;
4550          handler_->message() << CoinMessageEol;
4551     }
4552#if 0
4553     printf("IT %d %g %g(%d) %g(%d)\n",
4554            numberIterations_, objectiveValue(),
4555            sumPrimalInfeasibilities_, numberPrimalInfeasibilities_,
4556            sumDualInfeasibilities_, numberDualInfeasibilities_);
4557#endif
4558     double approximateObjective = objectiveValue_;
4559#ifdef CLP_REPORT_PROGRESS
4560     if (ixxxxxx >= ixxyyyy - 4 && ixxxxxx <= ixxyyyy) {
4561          char temp[20];
4562          sprintf(temp, "x_sol%d.out", ixxxxxx);
4563          FILE * fp = fopen(temp, "w");
4564          int nTotal = numberRows_ + numberColumns_;
4565          for (int i = 0; i < nTotal; i++)
4566               fprintf(fp, "%d %d %g %g %g %g %g\n",
4567                       i, status_[i], lower_[i], solution_[i], upper_[i], cost_[i], dj_[i]);
4568          fclose(fp);
4569          if (ixxxxxx == ixxyyyy)
4570               exit(6);
4571     }
4572#endif
4573     realDualInfeasibilities = sumDualInfeasibilities_;
4574     double saveTolerance = dualTolerance_;
4575     // If we need to carry on cleaning variables
4576     if (!numberPrimalInfeasibilities_ && (specialOptions_ & 1024) != 0 && CLEAN_FIXED) {
4577          for (int iRow = 0; iRow < numberRows_; iRow++) {
4578               int iPivot = pivotVariable_[iRow];
4579               if (!flagged(iPivot) && pivoted(iPivot)) {
4580                    // carry on
4581                    numberPrimalInfeasibilities_ = -1;
4582                    sumOfRelaxedPrimalInfeasibilities_ = 1.0;
4583                    sumPrimalInfeasibilities_ = 1.0;
4584                    break;
4585               }
4586          }
4587     }
4588     /* If we are primal feasible and any dual infeasibilities are on
4589        free variables then it is better to go to primal */
4590     if (!numberPrimalInfeasibilities_ && !numberDualInfeasibilitiesWithoutFree_ &&
4591               numberDualInfeasibilities_)
4592          problemStatus_ = 10;
4593     // dual bound coming in
4594     //double saveDualBound = dualBound_;
4595     bool needCleanFake = false;
4596     while (problemStatus_ <= -3) {
4597          int cleanDuals = 0;
4598          if (situationChanged != 0)
4599               cleanDuals = 1;
4600          int numberChangedBounds = 0;
4601          int doOriginalTolerance = 0;
4602          if ( lastCleaned == numberIterations_)
4603               doOriginalTolerance = 1;
4604          // check optimal
4605          // give code benefit of doubt
4606          if (sumOfRelaxedDualInfeasibilities_ == 0.0 &&
4607                    sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
4608               // say optimal (with these bounds etc)
4609               numberDualInfeasibilities_ = 0;
4610               sumDualInfeasibilities_ = 0.0;
4611               numberPrimalInfeasibilities_ = 0;
4612               sumPrimalInfeasibilities_ = 0.0;
4613          }
4614          //if (dualFeasible()||problemStatus_==-4||(primalFeasible()&&!numberDualInfeasibilitiesWithoutFree_)) {
4615          if (dualFeasible() || problemStatus_ == -4) {
4616               progress_.modifyObjective(objectiveValue_
4617                                         - bestPossibleImprovement_);
4618#ifdef COIN_DEVELOP
4619               if (sumDualInfeasibilities_ || bestPossibleImprovement_)
4620                    printf("improve %g dualinf %g -> %g\n",
4621                           bestPossibleImprovement_, sumDualInfeasibilities_,
4622                           sumDualInfeasibilities_ * dualBound_);
4623#endif
4624               // see if cutoff reached
4625               double limit = 0.0;
4626               getDblParam(ClpDualObjectiveLimit, limit);
4627#if 0
4628               if(fabs(limit) < 1.0e30 && objectiveValue()*optimizationDirection_ >
4629                         limit + 1.0e-7 + 1.0e-8 * fabs(limit) && !numberAtFakeBound()) {
4630                    //looks infeasible on objective
4631                    if (perturbation_ == 101) {
4632                         cleanDuals = 1;
4633                         // Save costs
4634                         int numberTotal = numberRows_ + numberColumns_;
4635                         double * saveCost = CoinCopyOfArray(cost_, numberTotal);
4636                         // make sure fake bounds are back
4637                         changeBounds(1, NULL, changeCost);
4638                         createRim4(false);
4639                         // make sure duals are current
4640                         computeDuals(givenDuals);
4641                         checkDualSolution();
4642                         if(objectiveValue()*optimizationDirection_ >
4643                                   limit + 1.0e-7 + 1.0e-8 * fabs(limit) && !numberDualInfeasibilities_) {
4644                              perturbation_ = 102; // stop any perturbations
4645                              printf("cutoff test succeeded\n");
4646                         } else {
4647                              printf("cutoff test failed\n");
4648                              // put back
4649                              memcpy(cost_, saveCost, numberTotal * sizeof(double));
4650                              // make sure duals are current
4651                              computeDuals(givenDuals);
4652                              checkDualSolution();
4653                              progress_.modifyObjective(-COIN_DBL_MAX);
4654                              problemStatus_ = -1;
4655                         }
4656                         delete [] saveCost;
4657                    }
4658               }
4659#endif
4660               if (primalFeasible() && !givenDuals) {
4661                    // may be optimal - or may be bounds are wrong
4662                    handler_->message(CLP_DUAL_BOUNDS, messages_)
4663                              << dualBound_
4664                              << CoinMessageEol;
4665                    // save solution in case unbounded
4666                    double * saveColumnSolution = NULL;
4667                    double * saveRowSolution = NULL;
4668                    bool inCbc = (specialOptions_ & (0x01000000 | 16384)) != 0;
4669                    if (!inCbc) {
4670                         saveColumnSolution = CoinCopyOfArray(columnActivityWork_, numberColumns_);
4671                         saveRowSolution = CoinCopyOfArray(rowActivityWork_, numberRows_);
4672                    }
4673                    numberChangedBounds = changeBounds(0, rowArray_[3], changeCost);
4674                    if (numberChangedBounds <= 0 && !numberDualInfeasibilities_) {
4675                         //looks optimal - do we need to reset tolerance
4676                         if (perturbation_ == 101) {
4677                              perturbation_ = 102; // stop any perturbations
4678                              cleanDuals = 1;
4679                              // make sure fake bounds are back
4680                              //computeObjectiveValue();
4681                              changeBounds(1, NULL, changeCost);
4682                              //computeObjectiveValue();
4683                              createRim4(false);
4684                              // make sure duals are current
4685                              computeDuals(givenDuals);
4686                              checkDualSolution();
4687                              progress_.modifyObjective(-COIN_DBL_MAX);
4688#define DUAL_TRY_FASTER
4689#ifdef DUAL_TRY_FASTER
4690                              if (numberDualInfeasibilities_) {
4691#endif
4692                                   numberChanged_ = 1; // force something to happen
4693                                   lastCleaned = numberIterations_ - 1;
4694#ifdef DUAL_TRY_FASTER
4695                              } else {
4696                                   //double value = objectiveValue_;
4697                                   computeObjectiveValue(true);
4698                                   //printf("old %g new %g\n",value,objectiveValue_);
4699                                   //numberChanged_=1;
4700                              }
4701#endif
4702                         }
4703                         if (lastCleaned < numberIterations_ && numberTimesOptimal_ < 4 &&
4704                                   (numberChanged_ || (specialOptions_ & 4096) == 0)) {
4705#if CLP_CAN_HAVE_ZERO_OBJ
4706                           if ((specialOptions_&2097152)==0) {
4707#endif
4708                              doOriginalTolerance = 2;
4709                              numberTimesOptimal_++;
4710                              changeMade_++; // say something changed
4711                              if (numberTimesOptimal_ == 1) {
4712                                   dualTolerance_ = dblParam_[ClpDualTolerance];
4713                              } else {
4714                                   if (numberTimesOptimal_ == 2) {
4715                                        // better to have small tolerance even if slower
4716                                        factorization_->zeroTolerance(CoinMin(factorization_->zeroTolerance(), 1.0e-15));
4717                                   }
4718                                   dualTolerance_ = dblParam_[ClpDualTolerance];
4719                                   dualTolerance_ *= pow(2.0, numberTimesOptimal_ - 1);
4720                              }
4721                              cleanDuals = 2; // If nothing changed optimal else primal
4722#if CLP_CAN_HAVE_ZERO_OBJ
4723                           } else {
4724                             // no cost - skip checks
4725                             problemStatus_=0;
4726                           }
4727#endif
4728                         } else {
4729                              problemStatus_ = 0; // optimal
4730                              if (lastCleaned < numberIterations_ && numberChanged_) {
4731                                   handler_->message(CLP_SIMPLEX_GIVINGUP, messages_)
4732                                             << CoinMessageEol;
4733                              }
4734                         }
4735                    } else {
4736                         cleanDuals = 1;
4737                         if (doOriginalTolerance == 1) {
4738                              // check unbounded
4739                              // find a variable with bad dj
4740                              int iSequence;
4741                              int iChosen = -1;
4742                              if (!inCbc) {
4743                                   double largest = 100.0 * primalTolerance_;
4744                                   for (iSequence = 0; iSequence < numberRows_ + numberColumns_;
4745                                             iSequence++) {
4746                                        double djValue = dj_[iSequence];
4747                                        double originalLo = originalLower(iSequence);
4748                                        double originalUp = originalUpper(iSequence);
4749                                        if (fabs(djValue) > fabs(largest)) {
4750                                             if (getStatus(iSequence) != basic) {
4751                                                  if (djValue > 0 && originalLo < -1.0e20) {
4752                                                       if (djValue > fabs(largest)) {
4753                                                            largest = djValue;
4754                                                            iChosen = iSequence;
4755                                                       }
4756                                                  } else if (djValue < 0 && originalUp > 1.0e20) {
4757                                                       if (-djValue > fabs(largest)) {
4758                                                            largest = djValue;
4759                                                            iChosen = iSequence;
4760                                                       }
4761                                                  }
4762                                             }
4763                                        }
4764                                   }
4765                              }
4766                              if (iChosen >= 0) {
4767                                   int iSave = sequenceIn_;
4768                                   sequenceIn_ = iChosen;
4769                                   unpack(rowArray_[1]);
4770                                   sequenceIn_ = iSave;
4771                                   // if dual infeasibilities then must be free vector so add in dual
4772                                   if (numberDualInfeasibilities_) {
4773                                        if (fabs(changeCost) > 1.0e-5)
4774                                          COIN_DETAIL_PRINT(printf("Odd free/unbounded combo\n"));
4775                                        changeCost += cost_[iChosen];
4776                                   }
4777                                   problemStatus_ = checkUnbounded(rowArray_[1], rowArray_[0],
4778                                                                   changeCost);
4779                                   rowArray_[1]->clear();
4780                              } else {
4781                                   problemStatus_ = -3;
4782                              }
4783                              if (problemStatus_ == 2 && perturbation_ == 101) {
4784                                   perturbation_ = 102; // stop any perturbations
4785                                   cleanDuals = 1;
4786                                   createRim4(false);
4787                                   progress_.modifyObjective(-COIN_DBL_MAX);
4788                                   problemStatus_ = -1;
4789                              }
4790                              if (problemStatus_ == 2) {
4791                                   // it is unbounded - restore solution
4792                                   // but first add in changes to non-basic
4793                                   int iColumn;
4794                                   double * original = columnArray_[0]->denseVector();
4795                                   for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
4796                                        if(getColumnStatus(iColumn) != basic)
4797                                             ray_[iColumn] +=
4798                                                  saveColumnSolution[iColumn] - original[iColumn];
4799                                        columnActivityWork_[iColumn] = original[iColumn];
4800                                   }
4801                                   CoinMemcpyN(saveRowSolution, numberRows_,
4802                                               rowActivityWork_);
4803                              }
4804                         } else {
4805                              doOriginalTolerance = 2;
4806                              rowArray_[0]->clear();
4807                         }
4808                    }
4809                    delete [] saveColumnSolution;
4810                    delete [] saveRowSolution;
4811               }
4812               if (problemStatus_ == -4 || problemStatus_ == -5) {
4813                    // may be infeasible - or may be bounds are wrong
4814                    numberChangedBounds = changeBounds(0, NULL, changeCost);
4815                    needCleanFake = true;
4816                    /* Should this be here as makes no difference to being feasible.
4817                       But seems to make a difference to run times. */
4818                    if (perturbation_ == 101 && 0) {
4819                         perturbation_ = 102; // stop any perturbations
4820                         cleanDuals = 1;
4821                         numberChangedBounds = 1;
4822                         // make sure fake bounds are back
4823                         changeBounds(1, NULL, changeCost);
4824                         needCleanFake = true;
4825                         createRim4(false);
4826                         progress_.modifyObjective(-COIN_DBL_MAX);
4827                    }
4828                    if ((numberChangedBounds <= 0 || dualBound_ > 1.0e20 ||
4829                              (largestPrimalError_ > 1.0 && dualBound_ > 1.0e17)) &&
4830                              (numberPivots < 4 || sumPrimalInfeasibilities_ > 1.0e-6)) {
4831                         problemStatus_ = 1; // infeasible
4832                         if (perturbation_ == 101) {
4833                              perturbation_ = 102; // stop any perturbations
4834                              //cleanDuals=1;
4835                              //numberChangedBounds=1;
4836                              //createRim4(false);
4837                         }
4838                         // but double check
4839                         if (!numberPrimalInfeasibilities_) {
4840                           problemStatus_=-1;
4841                           doOriginalTolerance=2;
4842                         }
4843                    } else {
4844                         problemStatus_ = -1; //iterate
4845                         cleanDuals = 1;
4846                         if (numberChangedBounds <= 0)
4847                              doOriginalTolerance = 2;
4848                         // and delete ray which has been created
4849                         delete [] ray_;
4850                         ray_ = NULL;
4851                    }
4852
4853               }
4854          } else {
4855               cleanDuals = 1;
4856          }
4857          if (problemStatus_ < 0) {
4858               if (doOriginalTolerance == 2) {
4859                    // put back original tolerance
4860                    lastCleaned = numberIterations_;
4861                    numberChanged_ = 0; // Number of variables with changed costs
4862                    handler_->message(CLP_DUAL_ORIGINAL, messages_)
4863                              << CoinMessageEol;
4864                    perturbation_ = 102; // stop any perturbations
4865#if 0
4866                    double * xcost = new double[numberRows_+numberColumns_];
4867                    double * xlower = new double[numberRows_+numberColumns_];
4868                    double * xupper = new double[numberRows_+numberColumns_];
4869                    double * xdj = new double[numberRows_+numberColumns_];
4870                    double * xsolution = new double[numberRows_+numberColumns_];
4871                    CoinMemcpyN(cost_, (numberRows_ + numberColumns_), xcost);
4872                    CoinMemcpyN(lower_, (numberRows_ + numberColumns_), xlower);
4873                    CoinMemcpyN(upper_, (numberRows_ + numberColumns_), xupper);
4874                    CoinMemcpyN(dj_, (numberRows_ + numberColumns_), xdj);
4875                    CoinMemcpyN(solution_, (numberRows_ + numberColumns_), xsolution);
4876#endif
4877                    createRim4(false);
4878                    progress_.modifyObjective(-COIN_DBL_MAX);
4879                    // make sure duals are current
4880                    computeDuals(givenDuals);
4881                    checkDualSolution();
4882#if 0
4883                    int i;
4884                    for (i = 0; i < numberRows_ + numberColumns_; i++) {
4885                         if (cost_[i] != xcost[i])
4886                              printf("** %d old cost %g new %g sol %g\n",
4887                                     i, xcost[i], cost_[i], solution_[i]);
4888                         if (lower_[i] != xlower[i])
4889                              printf("** %d old lower %g new %g sol %g\n",
4890                                     i, xlower[i], lower_[i], solution_[i]);
4891                         if (upper_[i] != xupper[i])
4892                              printf("** %d old upper %g new %g sol %g\n",
4893                                     i, xupper[i], upper_[i], solution_[i]);
4894                         if (dj_[i] != xdj[i])
4895                              printf("** %d old dj %g new %g sol %g\n",
4896                                     i, xdj[i], dj_[i], solution_[i]);
4897                         if (solution_[i] != xsolution[i])
4898                              printf("** %d old solution %g new %g sol %g\n",
4899                                     i, xsolution[i], solution_[i], solution_[i]);
4900                    }
4901                    //delete [] xcost;
4902                    //delete [] xupper;
4903                    //delete [] xlower;
4904                    //delete [] xdj;
4905                    //delete [] xsolution;
4906#endif
4907                    // put back bounds as they were if was optimal
4908                    if (doOriginalTolerance == 2 && cleanDuals != 2) {
4909                         changeMade_++; // say something changed
4910                         /* We may have already changed some bounds in this function
4911                            so save numberFake_ and add in.
4912
4913                            Worst that can happen is that we waste a bit of time  - but it must be finite.
4914                         */
4915                         //int saveNumberFake = numberFake_;
4916                         //resetFakeBounds(-1);
4917                         changeBounds(3, NULL, changeCost);
4918                         needCleanFake = true;
4919                         //numberFake_ += saveNumberFake;
4920                         //resetFakeBounds(-1);
4921                         cleanDuals = 2;
4922                         //cleanDuals=1;
4923                    }
4924#if 0
4925                    //int i;
4926                    for (i = 0; i < numberRows_ + numberColumns_; i++) {
4927                         if (cost_[i] != xcost[i])
4928                              printf("** %d old cost %g new %g sol %g\n",
4929                                     i, xcost[i], cost_[i], solution_[i]);
4930                         if (lower_[i] != xlower[i])
4931                              printf("** %d old lower %g new %g sol %g\n",
4932                                     i, xlower[i], lower_[i], solution_[i]);
4933                         if (upper_[i] != xupper[i])
4934                              printf("** %d old upper %g new %g sol %g\n",
4935                                     i, xupper[i], upper_[i], solution_[i]);
4936                         if (dj_[i] != xdj[i])
4937                              printf("** %d old dj %g new %g sol %g\n",
4938                                     i, xdj[i], dj_[i], solution_[i]);
4939                         if (solution_[i] != xsolution[i])
4940                              printf("** %d old solution %g new %g sol %g\n",
4941                                     i, xsolution[i], solution_[i], solution_[i]);
4942                    }
4943                    delete [] xcost;
4944                    delete [] xupper;
4945                    delete [] xlower;
4946                    delete [] xdj;
4947                    delete [] xsolution;
4948#endif
4949               }
4950               if (cleanDuals == 1 || (cleanDuals == 2 && !numberDualInfeasibilities_)) {
4951                    // make sure dual feasible
4952                    // look at all rows and columns
4953                    rowArray_[0]->clear();
4954                    columnArray_[0]->clear();
4955                    double objectiveChange = 0.0;
4956                    double savePrimalInfeasibilities = sumPrimalInfeasibilities_;
4957                    if (!numberIterations_) {
4958                      int nTotal = numberRows_ + numberColumns_;
4959                      if (arraysNotCreated) {
4960                        // create save arrays
4961                        delete [] saveStatus_;
4962                        delete [] savedSolution_;
4963                        saveStatus_ = new unsigned char [nTotal];
4964                        savedSolution_ = new double [nTotal];
4965                        arraysNotCreated = false;
4966                      }
4967                      // save arrays
4968                      CoinMemcpyN(status_, nTotal, saveStatus_);
4969                      CoinMemcpyN(rowActivityWork_,
4970                                  numberRows_, savedSolution_ + numberColumns_);
4971                      CoinMemcpyN(columnActivityWork_, numberColumns_, savedSolution_);
4972                    }
4973#if 0
4974                    double * xcost = new double[numberRows_+numberColumns_];
4975                    double * xlower = new double[numberRows_+numberColumns_];
4976                    double * xupper = new double[numberRows_+numberColumns_];
4977                    double * xdj = new double[numberRows_+numberColumns_];
4978                    double * xsolution = new double[numberRows_+numberColumns_];
4979                    CoinMemcpyN(cost_, (numberRows_ + numberColumns_), xcost);
4980                    CoinMemcpyN(lower_, (numberRows_ + numberColumns_), xlower);
4981                    CoinMemcpyN(upper_, (numberRows_ + numberColumns_), xupper);
4982                    CoinMemcpyN(dj_, (numberRows_ + numberColumns_), xdj);
4983                    CoinMemcpyN(solution_, (numberRows_ + numberColumns_), xsolution);
4984#endif
4985                    if (givenDuals)
4986                         dualTolerance_ = 1.0e50;
4987#if CLP_CAN_HAVE_ZERO_OBJ>1
4988                    if ((specialOptions_&2097152)==0) {
4989#endif
4990                    updateDualsInDual(rowArray_[0], columnArray_[0], rowArray_[1],
4991                                      0.0, objectiveChange, true);
4992#if CLP_CAN_HAVE_ZERO_OBJ>1
4993                    } else {
4994                      rowArray_[0]->clear();
4995                      rowArray_[1]->clear();
4996                      columnArray_[0]->clear();
4997                    }
4998#endif
4999                    dualTolerance_ = saveTolerance;
5000#if 0
5001                    int i;
5002                    for (i = 0; i < numberRows_ + numberColumns_; i++) {
5003                         if (cost_[i] != xcost[i])
5004                              printf("** %d old cost %g new %g sol %g\n",
5005                                     i, xcost[i], cost_[i], solution_[i]);
5006                         if (lower_[i] != xlower[i])
5007                              printf("** %d old lower %g new %g sol %g\n",
5008                                     i, xlower[i], lower_[i], solution_[i]);
5009                         if (upper_[i] != xupper[i])
5010                              printf("** %d old upper %g new %g sol %g\n",
5011                                     i, xupper[i], upper_[i], solution_[i]);
5012                         if (dj_[i] != xdj[i])
5013                              printf("** %d old dj %g new %g sol %g\n",
5014                                     i, xdj[i], dj_[i], solution_[i]);
5015                         if (solution_[i] != xsolution[i])
5016                              printf("** %d old solution %g new %g sol %g\n",
5017                                     i, xsolution[i], solution_[i], solution_[i]);
5018                    }
5019                    delete [] xcost;
5020                    delete [] xupper;
5021                    delete [] xlower;
5022                    delete [] xdj;
5023                    delete [] xsolution;
5024#endif
5025                    // for now - recompute all
5026                    gutsOfSolution(NULL, NULL);
5027                    if (givenDuals)
5028                         dualTolerance_ = 1.0e50;
5029#if CLP_CAN_HAVE_ZERO_OBJ>1
5030                    if ((specialOptions_&2097152)==0) {
5031#endif
5032                    updateDualsInDual(rowArray_[0], columnArray_[0], rowArray_[1],
5033                                      0.0, objectiveChange, true);
5034#if CLP_CAN_HAVE_ZERO_OBJ>1
5035                    } else {
5036                      rowArray_[0]->clear();
5037                      rowArray_[1]->clear();
5038                      columnArray_[0]->clear();
5039                    }
5040#endif
5041                    dualTolerance_ = saveTolerance;
5042                    if (!numberIterations_ && sumPrimalInfeasibilities_ >
5043                        1.0e5*(savePrimalInfeasibilities+1.0e3) &&
5044                        (moreSpecialOptions_ & (256|8192)) == 0) {
5045                      // Use primal
5046                      int nTotal = numberRows_ + numberColumns_;
5047                      CoinMemcpyN(saveStatus_, nTotal, status_);
5048                      CoinMemcpyN(savedSolution_ + numberColumns_ ,
5049                                  numberRows_, rowActivityWork_);
5050                      CoinMemcpyN(savedSolution_ ,
5051                                  numberColumns_, columnActivityWork_);
5052                      problemStatus_ = 10;
5053                      situationChanged = 0;
5054                    }
5055                    //assert(numberDualInfeasibilitiesWithoutFree_==0);
5056                    if (numberDualInfeasibilities_) {
5057                        if ((numberPrimalInfeasibilities_ || numberPivots)
5058                            && problemStatus_!=10) {
5059                              problemStatus_ = -1; // carry on as normal
5060                         } else {
5061                              problemStatus_ = 10; // try primal
5062#if COIN_DEVELOP>1
5063                              printf("returning at %d\n", __LINE__);
5064#endif
5065                         }
5066                    } else if (situationChanged == 2) {
5067                         problemStatus_ = -1; // carry on as normal
5068                         // need to reset bounds
5069                         changeBounds(3, NULL, changeCost);
5070                    }
5071                    situationChanged = 0;
5072               } else {
5073                    // iterate
5074                    if (cleanDuals != 2) {
5075                         problemStatus_ = -1;
5076                    } else {
5077                         problemStatus_ = 10; // try primal
5078#if COIN_DEVELOP>2
5079                         printf("returning at %d\n", __LINE__);
5080#endif
5081                    }
5082               }
5083          }
5084     }
5085     // unflag all variables (we may want to wait a bit?)
5086     if ((tentativeStatus != -2 && tentativeStatus != -1) && unflagVariables) {
5087          int iRow;
5088          int numberFlagged = 0;
5089          for (iRow = 0; iRow < numberRows_; iRow++) {
5090               int iPivot = pivotVariable_[iRow];
5091               if (flagged(iPivot)) {
5092                    numberFlagged++;
5093                    clearFlagged(iPivot);
5094               }
5095          }
5096#ifdef COIN_DEVELOP
5097          if (numberFlagged) {
5098               printf("unflagging %d variables - tentativeStatus %d probStat %d ninf %d nopt %d\n", numberFlagged, tentativeStatus,
5099                      problemStatus_, numberPrimalInfeasibilities_,
5100                      numberTimesOptimal_);
5101          }
5102#endif
5103          unflagVariables = numberFlagged > 0;
5104          if (numberFlagged && !numberPivots) {
5105               /* looks like trouble as we have not done any iterations.
5106               Try changing pivot tolerance then give it a few goes and give up */
5107               if (factorization_->pivotTolerance() < 0.9) {
5108                    factorization_->pivotTolerance(0.99);
5109                    problemStatus_ = -1;
5110               } else if (numberTimesOptimal_ < 3) {
5111                    numberTimesOptimal_++;
5112                    problemStatus_ = -1;
5113               } else {
5114                    unflagVariables = false;
5115                    //secondaryStatus_ = 1; // and say probably infeasible
5116                    if ((moreSpecialOptions_ & 256) == 0) {
5117                         // try primal
5118                         problemStatus_ = 10;
5119                    } else {
5120                         // almost certainly infeasible
5121                         problemStatus_ = 1;
5122                    }
5123#if COIN_DEVELOP>1
5124                    printf("returning at %d\n", __LINE__);
5125#endif
5126               }
5127          }
5128     }
5129     if (problemStatus_ < 0) {
5130          if (needCleanFake) {
5131               double dummyChangeCost = 0.0;
5132               changeBounds(3, NULL, dummyChangeCost);
5133          }
5134#if 0
5135          if (objectiveValue_ < lastObjectiveValue_ - 1.0e-8 *
5136                    CoinMax(fabs(objectivevalue_), fabs(lastObjectiveValue_))) {
5137          } else {
5138               lastObjectiveValue_ = objectiveValue_;
5139          }
5140#endif
5141          if (type == 0 || type == 1) {
5142               if (!type && arraysNotCreated) {
5143                    // create save arrays
5144                    delete [] saveStatus_;
5145                    delete [] savedSolution_;
5146                    saveStatus_ = new unsigned char [numberRows_+numberColumns_];
5147                    savedSolution_ = new double [numberRows_+numberColumns_];
5148               }
5149               // save arrays
5150               CoinMemcpyN(status_, numberColumns_ + numberRows_, saveStatus_);
5151               CoinMemcpyN(rowActivityWork_,
5152                           numberRows_, savedSolution_ + numberColumns_);
5153               CoinMemcpyN(columnActivityWork_, numberColumns_, savedSolution_);
5154               // save extra stuff
5155               int dummy;
5156               matrix_->generalExpanded(this, 5, dummy);
5157          }
5158          if (weightsSaved) {
5159               // restore weights (if saved) - also recompute infeasibility list
5160               if (!reallyBadProblems && (largestPrimalError_ < 100.0 || numberPivots > 10)) {
5161                    if (tentativeStatus > -3)
5162                         dualRowPivot_->saveWeights(this, (type < 2) ? 2 : 4);
5163                    else
5164                         dualRowPivot_->saveWeights(this, 3);
5165               } else {
5166                    // reset weights or scale back
5167                    dualRowPivot_->saveWeights(this, 6);
5168               }
5169          } else if (weightsSaved2 && numberPrimalInfeasibilities_) {
5170               dualRowPivot_->saveWeights(this, 3);
5171          }
5172     }
5173     // see if cutoff reached
5174     double limit = 0.0;
5175     getDblParam(ClpDualObjectiveLimit, limit);
5176#if 0
5177     if(fabs(limit) < 1.0e30 && objectiveValue()*optimizationDirection_ >
5178               limit + 100.0) {
5179          printf("lim %g obj %g %g - wo perturb %g sum dual %g\n",
5180                 limit, objectiveValue_, objectiveValue(), computeInternalObjectiveValue(), sumDualInfeasibilities_);
5181     }
5182#endif
5183     if(fabs(limit) < 1.0e30 && objectiveValue()*optimizationDirection_ >
5184               limit && !numberAtFakeBound()) {
5185          bool looksInfeasible = !numberDualInfeasibilities_;
5186          if (objectiveValue()*optimizationDirection_ > limit + fabs(0.1 * limit) + 1.0e2 * sumDualInfeasibilities_ + 1.0e4 &&
5187                    sumDualInfeasibilities_ < largestDualError_ && numberIterations_ > 0.5 * numberRows_ + 1000)
5188               looksInfeasible = true;
5189          if (looksInfeasible) {
5190               // Even if not perturbed internal costs may have changed
5191               // be careful
5192               if (true || numberIterations_) {
5193                    if(computeInternalObjectiveValue() > limit) {
5194                         problemStatus_ = 1;
5195                         secondaryStatus_ = 1; // and say was on cutoff
5196                    }
5197               } else {
5198                    problemStatus_ = 1;
5199                    secondaryStatus_ = 1; // and say was on cutoff
5200               }
5201          }
5202     }
5203     // If we are in trouble and in branch and bound give up
5204     if ((specialOptions_ & 1024) != 0) {
5205          int looksBad = 0;
5206          if (largestPrimalError_ * largestDualError_ > 1.0e2) {
5207               looksBad = 1;
5208          } else if (largestPrimalError_ > 1.0e-2
5209                     && objectiveValue_ > CoinMin(1.0e15, 1.0e3 * limit)) {
5210               looksBad = 2;
5211          }
5212          if (looksBad) {
5213               if (factorization_->pivotTolerance() < 0.9) {
5214                    // up tolerance
5215                    factorization_->pivotTolerance(CoinMin(factorization_->pivotTolerance() * 1.05 + 0.02, 0.91));
5216               } else if (numberIterations_ > 10000) {
5217                    if (handler_->logLevel() > 2)
5218                         printf("bad dual - saying infeasible %d\n", looksBad);
5219                    problemStatus_ = 1;
5220                    secondaryStatus_ = 1; // and say was on cutoff
5221               } else if (largestPrimalError_ > 1.0e5) {
5222                    {
5223                      //int iBigB = -1;
5224                         double bigB = 0.0;
5225                         //int iBigN = -1;
5226                         double bigN = 0.0;
5227                         for (int i = 0; i < numberRows_ + numberColumns_; i++) {
5228                              double value = fabs(solution_[i]);
5229                              if (getStatus(i) == basic) {
5230                                   if (value > bigB) {
5231                                        bigB = value;
5232                                        //iBigB = i;
5233                                   }
5234                              } else {
5235                                   if (value > bigN) {
5236                                        bigN = value;
5237                                        //iBigN = i;
5238                                   }
5239                              }
5240                         }
5241#ifdef CLP_INVESTIGATE
5242                         if (bigB > 1.0e8 || bigN > 1.0e8) {
5243                              if (handler_->logLevel() > 0)
5244                                   printf("it %d - basic %d %g, nonbasic %d %g\n",
5245                                          numberIterations_, iBigB, bigB, iBigN, bigN);
5246                         }
5247#endif
5248                    }
5249#if COIN_DEVELOP!=2
5250                    if (handler_->logLevel() > 2)
5251#endif
5252                         printf("bad dual - going to primal %d %g\n", looksBad, largestPrimalError_);
5253                    allSlackBasis(true);
5254                    problemStatus_ = 10;
5255               }
5256          }
5257     }
5258     if (problemStatus_ < 0 && !changeMade_) {