# source:stable/1.15/Clp/src/ClpSimplexDual.cpp@2006

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

changes for parallel and idiot

• Property svn:eol-style set to `native`
• Property svn:keywords set to `Id`
File size: 330.2 KB
Line
1/* \$Id: ClpSimplexDual.cpp 2006 2013-12-12 15:40:41Z forrest \$ */
2// Copyright (C) 2002, International Business Machines
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)
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_);
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_);
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                    // may have crept through - so may be optimal
1931                    // check any flagged variables
1932                    int iRow;
1933                    for (iRow = 0; iRow < numberRows_; iRow++) {
1934                         int iPivot = pivotVariable_[iRow];
1935                         if (flagged(iPivot))
1936                              break;
1937                    }
1938                    if (iRow < numberRows_ && numberPivots) {
1939                         // try factorization
1940                         returnCode = -2;
1941                    }
1942
1943                    if (useNumberFake || numberDualInfeasibilities_) {
1944                         // may be dual infeasible
1945                         if ((specialOptions_ & 1024) == 0)
1946                              problemStatus_ = -5;
1947                         else if (!useNumberFake && numberPrimalInfeasibilities_
1948                                   && !numberPivots)
1949                              problemStatus_ = 1;
1950                    } else {
1951                         if (iRow < numberRows_) {
1952#ifdef COIN_DEVELOP
1953                              std::cout << "Flagged variables at end - infeasible?" << std::endl;
1954                              printf("Probably infeasible - pivot was %g\n", alpha_);
1955#endif
1956                              //if (fabs(alpha_)<1.0e-4) {
1957                              //problemStatus_=1;
1958                              //} else {
1959#ifdef CLP_DEBUG
1960                              abort();
1961#endif
1962                              //}
1963                              problemStatus_ = -5;
1964                         } else {
1965                              problemStatus_ = 0;
1966#ifndef CLP_CHECK_NUMBER_PIVOTS
1967#define CLP_CHECK_NUMBER_PIVOTS 10
1968#endif
1969#if CLP_CHECK_NUMBER_PIVOTS < 20
1970                              if (numberPivots > CLP_CHECK_NUMBER_PIVOTS) {
1971#ifndef NDEBUG_CLP
1972                                   int nTotal = numberRows_ + numberColumns_;
1973                                   double * comp = CoinCopyOfArray(solution_, nTotal);
1974#endif
1975                                   computePrimals(rowActivityWork_, columnActivityWork_);
1976#ifndef NDEBUG_CLP
1977                                   double largest = 1.0e-5;
1978                                   int bad = -1;
1979                                   for (int i = 0; i < nTotal; i++) {
1980                                        double value = solution_[i];
1981                                        double larger = CoinMax(fabs(value), fabs(comp[i]));
1982                                        double tol = 1.0e-5 + 1.0e-5 * larger;
1983                                        double diff = fabs(value - comp[i]);
1984                                        if (diff - tol > largest) {
1985                                             bad = i;
1986                                             largest = diff - tol;
1987                                        }
1988                                   }
1989                                   if (bad >= 0)
1991#endif
1992                                   checkPrimalSolution(rowActivityWork_, columnActivityWork_);
1993                                   if (numberPrimalInfeasibilities_) {
1994#ifdef CLP_INVESTIGATE
1995                                        printf("XXX Infeas ? %d inf summing to %g\n", numberPrimalInfeasibilities_,
1996                                               sumPrimalInfeasibilities_);
1997#endif
1998                                        problemStatus_ = -1;
1999                                        returnCode = -2;
2000                                   }
2001#ifndef NDEBUG_CLP
2002                                   memcpy(solution_, comp, nTotal * sizeof(double));
2003                                   delete [] comp;
2004#endif
2005                              }
2006#endif
2007                              if (!problemStatus_) {
2008                                   // make it look OK
2009                                   numberPrimalInfeasibilities_ = 0;
2010                                   sumPrimalInfeasibilities_ = 0.0;
2011                                   numberDualInfeasibilities_ = 0;
2012                                   sumDualInfeasibilities_ = 0.0;
2013                                   // May be perturbed
2014                                   if (perturbation_ == 101 || numberChanged_) {
2015                                        numberChanged_ = 0; // Number of variables with changed costs
2016                                        perturbation_ = 102; // stop any perturbations
2017                                        //double changeCost;
2018                                        //changeBounds(1,NULL,changeCost);
2019                                        createRim4(false);
2020                                        // make sure duals are current
2021                                        computeDuals(givenDuals);
2022                                        checkDualSolution();
2023                                        progress_.modifyObjective(-COIN_DBL_MAX);
2024                                        if (numberDualInfeasibilities_) {
2025                                             problemStatus_ = 10; // was -3;
2026                                        } else {
2027                                             computeObjectiveValue(true);
2028                                        }
2029                                   } else if (numberPivots) {
2030                                        computeObjectiveValue(true);
2031                                   }
2032                                   if (numberPivots < -1000) {
2033                                        // objective may be wrong
2034                                        objectiveValue_ = innerProduct(cost_, numberColumns_ + numberRows_, solution_);
2035                                        objectiveValue_ += objective_->nonlinearOffset();
2036                                        objectiveValue_ /= (objectiveScale_ * rhsScale_);
2037                                        if ((specialOptions_ & 16384) == 0) {
2038                                             // and dual_ may be wrong (i.e. for fixed or basic)
2039                                             CoinIndexedVector * arrayVector = rowArray_[1];
2040                                             arrayVector->clear();
2041                                             int iRow;
2042                                             double * array = arrayVector->denseVector();
2043                                             /* Use dual_ instead of array
2044                                                Even though dual_ is only numberRows_ long this is
2045                                                okay as gets permuted to longer rowArray_[2]
2046                                             */
2047                                             arrayVector->setDenseVector(dual_);
2048                                             int * index = arrayVector->getIndices();
2049                                             int number = 0;
2050                                             for (iRow = 0; iRow < numberRows_; iRow++) {
2051                                                  int iPivot = pivotVariable_[iRow];
2052                                                  double value = cost_[iPivot];
2053                                                  dual_[iRow] = value;
2054                                                  if (value) {
2055                                                       index[number++] = iRow;
2056                                                  }
2057                                             }
2058                                             arrayVector->setNumElements(number);
2059                                             // Extended duals before "updateTranspose"
2060                                             matrix_->dualExpanded(this, arrayVector, NULL, 0);
2061                                             // Btran basic costs
2062                                             rowArray_[2]->clear();
2063                                             factorization_->updateColumnTranspose(rowArray_[2], arrayVector);
2064                                             // and return vector
2065                                             arrayVector->setDenseVector(array);
2066                                        }
2067                                   }
2068                                   sumPrimalInfeasibilities_ = 0.0;
2069                              }
2070                              if ((specialOptions_&(1024 + 16384)) != 0 && !problemStatus_) {
2071                                   CoinIndexedVector * arrayVector = rowArray_[1];
2072                                   arrayVector->clear();
2073                                   double * rhs = arrayVector->denseVector();
2074                                   times(1.0, solution_, rhs);
2075#ifdef CHECK_ACCURACY
2076                                   bool bad = false;
2077#endif
2078                                   bool bad2 = false;
2079                                   int i;
2080                                   for ( i = 0; i < numberRows_; i++) {
2081                                        if (rhs[i] < rowLowerWork_[i] - primalTolerance_ ||
2082                                                  rhs[i] > rowUpperWork_[i] + primalTolerance_) {
2083                                             bad2 = true;
2084#ifdef CHECK_ACCURACY
2085                                             printf("row %d out of bounds %g, %g correct %g bad %g\n", i,
2086                                                    rowLowerWork_[i], rowUpperWork_[i],
2087                                                    rhs[i], rowActivityWork_[i]);
2088#endif
2089                                        } else if (fabs(rhs[i] - rowActivityWork_[i]) > 1.0e-3) {
2090#ifdef CHECK_ACCURACY
2091                                             bad = true;
2092                                             printf("row %d correct %g bad %g\n", i, rhs[i], rowActivityWork_[i]);
2093#endif
2094                                        }
2095                                        rhs[i] = 0.0;
2096                                   }
2097                                   for ( i = 0; i < numberColumns_; i++) {
2098                                        if (solution_[i] < columnLowerWork_[i] - primalTolerance_ ||
2099                                                  solution_[i] > columnUpperWork_[i] + primalTolerance_) {
2100                                             bad2 = true;
2101#ifdef CHECK_ACCURACY
2102                                             printf("column %d out of bounds %g, %g correct %g bad %g\n", i,
2103                                                    columnLowerWork_[i], columnUpperWork_[i],
2104                                                    solution_[i], columnActivityWork_[i]);
2105#endif
2106                                        }
2107                                   }
2108                                   if (bad2) {
2109                                        problemStatus_ = -3;
2110                                        returnCode = -2;
2111                                        // Force to re-factorize early next time
2112                                        int numberPivots = factorization_->pivots();
2113                                        forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
2114                                   }
2115                              }
2116                         }
2117                    }
2118               } else {
2119                    problemStatus_ = -3;
2120                    returnCode = -2;
2121                    // Force to re-factorize early next time
2122                    int numberPivots = factorization_->pivots();
2123                    forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
2124               }
2125               break;
2126          }
2127     }
2128     if (givenDuals) {
2129          CoinMemcpyN(dj_, numberRows_ + numberColumns_, givenDuals);
2130          // get rid of any values pass array
2131          delete [] candidateList;
2132     }
2133     delete [] dubiousWeights;
2134#ifdef CLP_REPORT_PROGRESS
2135     if (ixxxxxx > ixxyyyy - 5) {
2136          int nTotal = numberColumns_ + numberRows_;
2137          double oldObj = 0.0;
2138          double newObj = 0.0;
2139          for (int i = 0; i < nTotal; i++) {
2140               if (savePSol[i])
2141                    oldObj += savePSol[i] * saveCost[i];
2142               if (solution_[i])
2143                    newObj += solution_[i] * cost_[i];
2144               bool printIt = false;
2145               if (cost_[i] != saveCost[i])
2146                    printIt = true;
2147               if (status_[i] != saveStat[i])
2148                    printIt = true;
2149               if (printIt)
2150                    printf("%d old %d cost %g sol %g, new %d cost %g sol %g\n",
2151                           i, saveStat[i], saveCost[i], savePSol[i],
2152                           status_[i], cost_[i], solution_[i]);
2153               // difference
2154               savePSol[i] = solution_[i] - savePSol[i];
2155          }
2156          printf("exit pivots %d, old obj %g new %g\n",
2157                 factorization_->pivots(),
2158                 oldObj, newObj);
2159          memset(saveDj, 0, numberRows_ * sizeof(double));
2160          times(1.0, savePSol, saveDj);
2161          double largest = 1.0e-6;
2162          int k = -1;
2163          for (int i = 0; i < numberRows_; i++) {
2164               saveDj[i] -= savePSol[i+numberColumns_];
2165               if (fabs(saveDj[i]) > largest) {
2166                    largest = fabs(saveDj[i]);
2167                    k = i;
2168               }
2169          }
2170          if (k >= 0)
2171               printf("Not null %d %g\n", k, largest);
2172     }
2173     delete [] savePSol ;
2174     delete [] saveDj ;
2175     delete [] saveCost ;
2176     delete [] saveStat ;
2177#endif
2178     return returnCode;
2179}
2180/* The duals are updated by the given arrays.
2181   Returns number of infeasibilities.
2182   rowArray and columnarray will have flipped
2183   The output vector has movement (row length array) */
2184int
2185ClpSimplexDual::updateDualsInDual(CoinIndexedVector * rowArray,
2186                                  CoinIndexedVector * columnArray,
2187                                  CoinIndexedVector * outputArray,
2188                                  double theta,
2189                                  double & objectiveChange,
2190                                  bool fullRecompute)
2191{
2192
2193     outputArray->clear();
2194
2195
2196     int numberInfeasibilities = 0;
2197     int numberRowInfeasibilities = 0;
2198
2199     // get a tolerance
2200     double tolerance = dualTolerance_;
2201     // we can't really trust infeasibilities if there is dual error
2202     double error = CoinMin(1.0e-2, largestDualError_);
2203     // allow tolerance at least slightly bigger than standard
2204     tolerance = tolerance +  error;
2205
2206     double changeObj = 0.0;
2207
2208     // Coding is very similar but we can save a bit by splitting
2209     // Do rows
2210     if (!fullRecompute) {
2211          int i;
2212          double * COIN_RESTRICT reducedCost = djRegion(0);
2213          const double * COIN_RESTRICT lower = lowerRegion(0);
2214          const double * COIN_RESTRICT upper = upperRegion(0);
2215          const double * COIN_RESTRICT cost = costRegion(0);
2216          double * COIN_RESTRICT work;
2217          int number;
2218          int * COIN_RESTRICT which;
2219          const unsigned char * COIN_RESTRICT statusArray = status_ + numberColumns_;
2220          assert(rowArray->packedMode());
2221          work = rowArray->denseVector();
2222          number = rowArray->getNumElements();
2223          which = rowArray->getIndices();
2224          double multiplier[] = { -1.0, 1.0};
2225          for (i = 0; i < number; i++) {
2226               int iSequence = which[i];
2227               double alphaI = work[i];
2228               work[i] = 0.0;
2229               int iStatus = (statusArray[iSequence] & 3) - 1;
2230               if (iStatus) {
2231                    double value = reducedCost[iSequence] - theta * alphaI;
2232                    reducedCost[iSequence] = value;
2233                    double mult = multiplier[iStatus-1];
2234                    value *= mult;
2235                    if (value < -tolerance) {
2236                         // flipping bounds
2237                         double movement = mult * (lower[iSequence] - upper[iSequence]);
2238                         which[numberInfeasibilities++] = iSequence;
2239#ifndef NDEBUG
2240                         if (fabs(movement) >= 1.0e30)
2241                              resetFakeBounds(-1000 - iSequence);
2242#endif
2243#ifdef CLP_DEBUG
2244                         if ((handler_->logLevel() & 32))
2245                              printf("%d %d, new dj %g, alpha %g, movement %g\n",
2246                                     0, iSequence, value, alphaI, movement);
2247#endif
2248                         changeObj -= movement * cost[iSequence];
2250                    }
2251               }
2252          }
2253          // Do columns
2254          reducedCost = djRegion(1);
2255          lower = lowerRegion(1);
2256          upper = upperRegion(1);
2257          cost = costRegion(1);
2258          // set number of infeasibilities in row array
2259          numberRowInfeasibilities = numberInfeasibilities;
2260          rowArray->setNumElements(numberInfeasibilities);
2261          numberInfeasibilities = 0;
2262          work = columnArray->denseVector();
2263          number = columnArray->getNumElements();
2264          which = columnArray->getIndices();
2265          if ((moreSpecialOptions_ & 8) != 0) {
2266               const unsigned char * COIN_RESTRICT statusArray = status_;
2267               for (i = 0; i < number; i++) {
2268                    int iSequence = which[i];
2269                    double alphaI = work[i];
2270                    work[i] = 0.0;
2271
2272                    int iStatus = (statusArray[iSequence] & 3) - 1;
2273                    if (iStatus) {
2274                         double value = reducedCost[iSequence] - theta * alphaI;
2275                         reducedCost[iSequence] = value;
2276                         double mult = multiplier[iStatus-1];
2277                         value *= mult;
2278                         if (value < -tolerance) {
2279                              // flipping bounds
2280                              double movement = mult * (upper[iSequence] - lower[iSequence]);
2281                              which[numberInfeasibilities++] = iSequence;
2282#ifndef NDEBUG
2283                              if (fabs(movement) >= 1.0e30)
2284                                   resetFakeBounds(-1000 - iSequence);
2285#endif
2286#ifdef CLP_DEBUG
2287                              if ((handler_->logLevel() & 32))
2288                                   printf("%d %d, new dj %g, alpha %g, movement %g\n",
2289                                          1, iSequence, value, alphaI, movement);
2290#endif
2291                              changeObj += movement * cost[iSequence];
2292                              matrix_->add(this, outputArray, iSequence, movement);
2293                         }
2294                    }
2295               }
2296          } else {
2297               for (i = 0; i < number; i++) {
2298                    int iSequence = which[i];
2299                    double alphaI = work[i];
2300                    work[i] = 0.0;
2301
2302                    Status status = getStatus(iSequence);
2303                    if (status == atLowerBound) {
2304                         double value = reducedCost[iSequence] - theta * alphaI;
2305                         reducedCost[iSequence] = value;
2306                         double movement = 0.0;
2307
2308                         if (value < -tolerance) {
2309                              // to upper bound
2310                              which[numberInfeasibilities++] = iSequence;
2311                              movement = upper[iSequence] - lower[iSequence];
2312#ifndef NDEBUG
2313                              if (fabs(movement) >= 1.0e30)
2314                                   resetFakeBounds(-1000 - iSequence);
2315#endif
2316#ifdef CLP_DEBUG
2317                              if ((handler_->logLevel() & 32))
2318                                   printf("%d %d, new dj %g, alpha %g, movement %g\n",
2319                                          1, iSequence, value, alphaI, movement);
2320#endif
2321                              changeObj += movement * cost[iSequence];
2322                              matrix_->add(this, outputArray, iSequence, movement);
2323                         }
2324                    } else if (status == atUpperBound) {
2325                         double value = reducedCost[iSequence] - theta * alphaI;
2326                         reducedCost[iSequence] = value;
2327                         double movement = 0.0;
2328
2329                         if (value > tolerance) {
2330                              // to lower bound (if swap)
2331                              which[numberInfeasibilities++] = iSequence;
2332                              movement = lower[iSequence] - upper[iSequence];
2333#ifndef NDEBUG
2334                              if (fabs(movement) >= 1.0e30)
2335                                   resetFakeBounds(-1000 - iSequence);
2336#endif
2337#ifdef CLP_DEBUG
2338                              if ((handler_->logLevel() & 32))
2339                                   printf("%d %d, new dj %g, alpha %g, movement %g\n",
2340                                          1, iSequence, value, alphaI, movement);
2341#endif
2342                              changeObj += movement * cost[iSequence];
2343                              matrix_->add(this, outputArray, iSequence, movement);
2344                         }
2345                    } else if (status == isFree) {
2346                         double value = reducedCost[iSequence] - theta * alphaI;
2347                         reducedCost[iSequence] = value;
2348                    }
2349               }
2350          }
2351     } else  {
2352          double * COIN_RESTRICT solution = solutionRegion(0);
2353          double * COIN_RESTRICT reducedCost = djRegion(0);
2354          const double * COIN_RESTRICT lower = lowerRegion(0);
2355          const double * COIN_RESTRICT upper = upperRegion(0);
2356          const double * COIN_RESTRICT cost = costRegion(0);
2357          int * COIN_RESTRICT which;
2358          which = rowArray->getIndices();
2359          int iSequence;
2360          for (iSequence = 0; iSequence < numberRows_; iSequence++) {
2361               double value = reducedCost[iSequence];
2362
2363               Status status = getStatus(iSequence + numberColumns_);
2364               // more likely to be at upper bound ?
2365               if (status == atUpperBound) {
2366                    double movement = 0.0;
2367                    //#define NO_SWAP7
2368                    if (value > tolerance) {
2369                         // to lower bound (if swap)
2370                         // put back alpha
2371                         which[numberInfeasibilities++] = iSequence;
2372                         movement = lower[iSequence] - upper[iSequence];
2373                         changeObj += movement * cost[iSequence];
2375#ifndef NDEBUG
2376                         if (fabs(movement) >= 1.0e30)
2377                              resetFakeBounds(-1000 - iSequence);
2378#endif
2379#ifndef NO_SWAP7
2380                    } else if (value > -tolerance) {
2381                         // at correct bound but may swap
2382                         FakeBound bound = getFakeBound(iSequence + numberColumns_);
2383                         if (bound == ClpSimplexDual::upperFake) {
2384                              movement = lower[iSequence] - upper[iSequence];
2385#ifndef NDEBUG
2386                              if (fabs(movement) >= 1.0e30)
2387                                   resetFakeBounds(-1000 - iSequence);
2388#endif
2389                              setStatus(iSequence + numberColumns_, atLowerBound);
2390                              solution[iSequence] = lower[iSequence];
2391                              changeObj += movement * cost[iSequence];
2392                              //numberFake_--;
2393                              //setFakeBound(iSequence+numberColumns_,noFake);
2394                         }
2395#endif
2396                    }
2397               } else if (status == atLowerBound) {
2398                    double movement = 0.0;
2399
2400                    if (value < -tolerance) {
2401                         // to upper bound
2402                         // put back alpha
2403                         which[numberInfeasibilities++] = iSequence;
2404                         movement = upper[iSequence] - lower[iSequence];
2405#ifndef NDEBUG
2406                         if (fabs(movement) >= 1.0e30)
2407                              resetFakeBounds(-1000 - iSequence);
2408#endif
2409                         changeObj += movement * cost[iSequence];
2411#ifndef NO_SWAP7
2412                    } else if (value < tolerance) {
2413                         // at correct bound but may swap
2414                         FakeBound bound = getFakeBound(iSequence + numberColumns_);
2415                         if (bound == ClpSimplexDual::lowerFake) {
2416                              movement = upper[iSequence] - lower[iSequence];
2417#ifndef NDEBUG
2418                              if (fabs(movement) >= 1.0e30)
2419                                   resetFakeBounds(-1000 - iSequence);
2420#endif
2421                              setStatus(iSequence + numberColumns_, atUpperBound);
2422                              solution[iSequence] = upper[iSequence];
2423                              changeObj += movement * cost[iSequence];
2424                              //numberFake_--;
2425                              //setFakeBound(iSequence+numberColumns_,noFake);
2426                         }
2427#endif
2428                    }
2429               }
2430          }
2431          // Do columns
2432          solution = solutionRegion(1);
2433          reducedCost = djRegion(1);
2434          lower = lowerRegion(1);
2435          upper = upperRegion(1);
2436          cost = costRegion(1);
2437          // set number of infeasibilities in row array
2438          numberRowInfeasibilities = numberInfeasibilities;
2439          rowArray->setNumElements(numberInfeasibilities);
2440          numberInfeasibilities = 0;
2441          which = columnArray->getIndices();
2442          for (iSequence = 0; iSequence < numberColumns_; iSequence++) {
2443               double value = reducedCost[iSequence];
2444
2445               Status status = getStatus(iSequence);
2446               if (status == atLowerBound) {
2447                    double movement = 0.0;
2448
2449                    if (value < -tolerance) {
2450                         // to upper bound
2451                         // put back alpha
2452                         which[numberInfeasibilities++] = iSequence;
2453                         movement = upper[iSequence] - lower[iSequence];
2454#ifndef NDEBUG
2455                         if (fabs(movement) >= 1.0e30)
2456                              resetFakeBounds(-1000 - iSequence);
2457#endif
2458                         changeObj += movement * cost[iSequence];
2459                         matrix_->add(this, outputArray, iSequence, movement);
2460#ifndef NO_SWAP7
2461                    } else if (value < tolerance) {
2462                         // at correct bound but may swap
2463                         FakeBound bound = getFakeBound(iSequence);
2464                         if (bound == ClpSimplexDual::lowerFake) {
2465                              movement = upper[iSequence] - lower[iSequence];
2466#ifndef NDEBUG
2467                              if (fabs(movement) >= 1.0e30)
2468                                   resetFakeBounds(-1000 - iSequence);
2469#endif
2470                              setStatus(iSequence, atUpperBound);
2471                              solution[iSequence] = upper[iSequence];
2472                              changeObj += movement * cost[iSequence];
2473                              //numberFake_--;
2474                              //setFakeBound(iSequence,noFake);
2475                         }
2476#endif
2477                    }
2478               } else if (status == atUpperBound) {
2479                    double movement = 0.0;
2480
2481                    if (value > tolerance) {
2482                         // to lower bound (if swap)
2483                         // put back alpha
2484                         which[numberInfeasibilities++] = iSequence;
2485                         movement = lower[iSequence] - upper[iSequence];
2486#ifndef NDEBUG
2487                         if (fabs(movement) >= 1.0e30)
2488                              resetFakeBounds(-1000 - iSequence);
2489#endif
2490                         changeObj += movement * cost[iSequence];
2491                         matrix_->add(this, outputArray, iSequence, movement);
2492#ifndef NO_SWAP7
2493                    } else if (value > -tolerance) {
2494                         // at correct bound but may swap
2495                         FakeBound bound = getFakeBound(iSequence);
2496                         if (bound == ClpSimplexDual::upperFake) {
2497                              movement = lower[iSequence] - upper[iSequence];
2498#ifndef NDEBUG
2499                              if (fabs(movement) >= 1.0e30)
2500                                   resetFakeBounds(-1000 - iSequence);
2501#endif
2502                              setStatus(iSequence, atLowerBound);
2503                              solution[iSequence] = lower[iSequence];
2504                              changeObj += movement * cost[iSequence];
2505                              //numberFake_--;
2506                              //setFakeBound(iSequence,noFake);
2507                         }
2508#endif
2509                    }
2510               }
2511          }
2512     }
2513
2514#ifdef CLP_DEBUG
2515     if (fullRecompute && numberFake_ && (handler_->logLevel() & 16) != 0)
2516          printf("%d fake after full update\n", numberFake_);
2517#endif
2518     // set number of infeasibilities
2519     columnArray->setNumElements(numberInfeasibilities);
2520     numberInfeasibilities += numberRowInfeasibilities;
2521     if (fullRecompute) {
2522          // do actual flips
2523          flipBounds(rowArray, columnArray);
2524     }
2525     objectiveChange += changeObj;
2526     return numberInfeasibilities;
2527}
2528void
2529ClpSimplexDual::updateDualsInValuesPass(CoinIndexedVector * rowArray,
2530                                        CoinIndexedVector * columnArray,
2531                                        double theta)
2532{
2533
2534     // use a tighter tolerance except for all being okay
2535     double tolerance = dualTolerance_;
2536
2537     // Coding is very similar but we can save a bit by splitting
2538     // Do rows
2539     {
2540          int i;
2541          double * reducedCost = djRegion(0);
2542          double * work;
2543          int number;
2544          int * which;
2545          work = rowArray->denseVector();
2546          number = rowArray->getNumElements();
2547          which = rowArray->getIndices();
2548          for (i = 0; i < number; i++) {
2549               int iSequence = which[i];
2550               double alphaI = work[i];
2551               double value = reducedCost[iSequence] - theta * alphaI;
2552               work[i] = 0.0;
2553               reducedCost[iSequence] = value;
2554
2555               Status status = getStatus(iSequence + numberColumns_);
2556               // more likely to be at upper bound ?
2557               if (status == atUpperBound) {
2558
2559                    if (value > tolerance)
2560                         reducedCost[iSequence] = 0.0;
2561               } else if (status == atLowerBound) {
2562
2563                    if (value < -tolerance) {
2564                         reducedCost[iSequence] = 0.0;
2565                    }
2566               }
2567          }
2568     }
2569     rowArray->setNumElements(0);
2570
2571     // Do columns
2572     {
2573          int i;
2574          double * reducedCost = djRegion(1);
2575          double * work;
2576          int number;
2577          int * which;
2578          work = columnArray->denseVector();
2579          number = columnArray->getNumElements();
2580          which = columnArray->getIndices();
2581
2582          for (i = 0; i < number; i++) {
2583               int iSequence = which[i];
2584               double alphaI = work[i];
2585               double value = reducedCost[iSequence] - theta * alphaI;
2586               work[i] = 0.0;
2587               reducedCost[iSequence] = value;
2588
2589               Status status = getStatus(iSequence);
2590               if (status == atLowerBound) {
2591                    if (value < -tolerance)
2592                         reducedCost[iSequence] = 0.0;
2593               } else if (status == atUpperBound) {
2594                    if (value > tolerance)
2595                         reducedCost[iSequence] = 0.0;
2596               }
2597          }
2598     }
2599     columnArray->setNumElements(0);
2600}
2601/*
2602   Chooses dual pivot row
2603   Would be faster with separate region to scan
2604   and will have this (with square of infeasibility) when steepest
2605   For easy problems we can just choose one of the first rows we look at
2606*/
2607void
2609{
2610     // get pivot row using whichever method it is
2611     int chosenRow = -1;
2612#ifdef FORCE_FOLLOW
2613     bool forceThis = false;
2614     if (!fpFollow && strlen(forceFile)) {
2615          fpFollow = fopen(forceFile, "r");
2616          assert (fpFollow);
2617     }
2618     if (fpFollow) {
2619          if (numberIterations_ <= force_iteration) {
2620               // read to next Clp0102
2621               char temp[300];
2622               while (fgets(temp, 250, fpFollow)) {
2623                    if (strncmp(temp, "Clp0102", 7))
2624                         continue;
2625                    char cin, cout;
2626                    sscanf(temp + 9, "%d%*f%*s%*c%c%d%*s%*c%c%d",
2627                           &force_iteration, &cin, &force_in, &cout, &force_out);
2628                    if (cin == 'R')
2629                         force_in += numberColumns_;
2630                    if (cout == 'R')
2631                         force_out += numberColumns_;
2632                    forceThis = true;
2633                    assert (numberIterations_ == force_iteration - 1);
2634                    printf("Iteration %d will force %d out and %d in\n",
2635                           force_iteration, force_out, force_in);
2636                    alreadyChosen = force_out;
2637                    break;
2638               }
2639          } else {
2640               // use old
2641               forceThis = true;
2642          }
2643          if (!forceThis) {
2644               fclose(fpFollow);
2645               fpFollow = NULL;
2646               forceFile = "";
2647          }
2648     }
2649#endif
2650     //double freeAlpha = 0.0;
2651     if (alreadyChosen < 0) {
2652          // first see if any free variables and put them in basis
2653          int nextFree = nextSuperBasic();
2654          //nextFree=-1; //off
2655          if (nextFree >= 0) {
2656               // unpack vector and find a good pivot
2657               unpack(rowArray_[1], nextFree);
2658               factorization_->updateColumn(rowArray_[2], rowArray_[1]);
2659
2660               double * work = rowArray_[1]->denseVector();
2661               int number = rowArray_[1]->getNumElements();
2662               int * which = rowArray_[1]->getIndices();
2663               double bestFeasibleAlpha = 0.0;
2664               int bestFeasibleRow = -1;
2665               double bestInfeasibleAlpha = 0.0;
2666               int bestInfeasibleRow = -1;
2667               int i;
2668
2669               for (i = 0; i < number; i++) {
2670                    int iRow = which[i];
2671                    double alpha = fabs(work[iRow]);
2672                    if (alpha > 1.0e-3) {
2673                         int iSequence = pivotVariable_[iRow];
2674                         double value = solution_[iSequence];
2675                         double lower = lower_[iSequence];
2676                         double upper = upper_[iSequence];
2677                         double infeasibility = 0.0;
2678                         if (value > upper)
2679                              infeasibility = value - upper;
2680                         else if (value < lower)
2681                              infeasibility = lower - value;
2682                         if (infeasibility * alpha > bestInfeasibleAlpha && alpha > 1.0e-1) {
2683                              if (!flagged(iSequence)) {
2684                                   bestInfeasibleAlpha = infeasibility * alpha;
2685                                   bestInfeasibleRow = iRow;
2686                              }
2687                         }
2688                         if (alpha > bestFeasibleAlpha && (lower > -1.0e20 || upper < 1.0e20)) {
2689                              bestFeasibleAlpha = alpha;
2690                              bestFeasibleRow = iRow;
2691                         }
2692                    }
2693               }
2694               if (bestInfeasibleRow >= 0)
2695                    chosenRow = bestInfeasibleRow;
2696               else if (bestFeasibleAlpha > 1.0e-2)
2697                    chosenRow = bestFeasibleRow;
2698               if (chosenRow >= 0) {
2699                    pivotRow_ = chosenRow;
2700                    //freeAlpha = work[chosenRow];
2701               }
2702               rowArray_[1]->clear();
2703          }
2704     } else {
2705          // in values pass
2706          chosenRow = alreadyChosen;
2707#ifdef FORCE_FOLLOW
2708          if(forceThis) {
2709               alreadyChosen = -1;
2710               chosenRow = -1;
2711               for (int i = 0; i < numberRows_; i++) {
2712                    if (pivotVariable_[i] == force_out) {
2713                         chosenRow = i;
2714                         break;
2715                    }
2716               }
2717               assert (chosenRow >= 0);
2718          }
2719#endif
2720          pivotRow_ = chosenRow;
2721     }
2722     if (chosenRow < 0)
2723          pivotRow_ = dualRowPivot_->pivotRow();
2724
2725     if (pivotRow_ >= 0) {
2726          sequenceOut_ = pivotVariable_[pivotRow_];
2727          valueOut_ = solution_[sequenceOut_];
2728          lowerOut_ = lower_[sequenceOut_];
2729          upperOut_ = upper_[sequenceOut_];
2730          if (alreadyChosen < 0) {
2731               // if we have problems we could try other way and hope we get a
2732               // zero pivot?
2733               if (valueOut_ > upperOut_) {
2734                    directionOut_ = -1;
2735                    dualOut_ = valueOut_ - upperOut_;
2736               } else if (valueOut_ < lowerOut_) {
2737                    directionOut_ = 1;
2738                    dualOut_ = lowerOut_ - valueOut_;
2739               } else {
2740#if 1
2741                    // odd (could be free) - it's feasible - go to nearest
2742                    if (valueOut_ - lowerOut_ < upperOut_ - valueOut_) {
2743                         directionOut_ = 1;
2744                         dualOut_ = lowerOut_ - valueOut_;
2745                    } else {
2746                         directionOut_ = -1;
2747                         dualOut_ = valueOut_ - upperOut_;
2748                    }
2749#else
2750                    // odd (could be free) - it's feasible - improve obj
2751                    printf("direction from alpha of %g is %d\n",
2752                           freeAlpha, freeAlpha > 0.0 ? 1 : -1);
2753                    if (valueOut_ - lowerOut_ > 1.0e20)
2754                         freeAlpha = 1.0;
2755                    else if(upperOut_ - valueOut_ > 1.0e20)
2756                         freeAlpha = -1.0;
2757                    //if (valueOut_-lowerOut_<upperOut_-valueOut_) {
2758                    if (freeAlpha < 0.0) {
2759                         directionOut_ = 1;
2760                         dualOut_ = lowerOut_ - valueOut_;
2761                    } else {
2762                         directionOut_ = -1;
2763                         dualOut_ = valueOut_ - upperOut_;
2764                    }
2765                    printf("direction taken %d - bounds %g %g %g\n",
2766                           directionOut_, lowerOut_, valueOut_, upperOut_);
2767#endif
2768               }
2769#ifdef CLP_DEBUG
2770               assert(dualOut_ >= 0.0);
2771#endif
2772          } else {
2773               // in values pass so just use sign of dj
2774               // We don't want to go through any barriers so set dualOut low
2775               // free variables will never be here
2776               dualOut_ = 1.0e-6;
2777               if (dj_[sequenceOut_] > 0.0) {
2778                    // this will give a -1 in pivot row (as slacks are -1.0)
2779                    directionOut_ = 1;
2780               } else {
2781                    directionOut_ = -1;
2782               }
2783          }
2784     }
2785     return ;
2786}
2787// Checks if any fake bounds active - if so returns number and modifies
2788// dualBound_ and everything.
2789// Free variables will be left as free
2790// Returns number of bounds changed if >=0
2791// Returns -1 if not initialize and no effect
2792// Fills in changeVector which can be used to see if unbounded
2793// and cost of change vector
2794int
2795ClpSimplexDual::changeBounds(int initialize,
2796                             CoinIndexedVector * outputArray,
2797                             double & changeCost)
2798{
2799     numberFake_ = 0;
2800     if (!initialize) {
2801          int numberInfeasibilities;
2802          double newBound;
2803          newBound = 5.0 * dualBound_;
2804          numberInfeasibilities = 0;
2805          changeCost = 0.0;
2806          // put back original bounds and then check
2807          createRim1(false);
2808          int iSequence;
2809          // bounds will get bigger - just look at ones at bounds
2810          for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
2811               double lowerValue = lower_[iSequence];
2812               double upperValue = upper_[iSequence];
2813               double value = solution_[iSequence];
2814               setFakeBound(iSequence, ClpSimplexDual::noFake);
2815               switch(getStatus(iSequence)) {
2816
2817               case basic:
2818               case ClpSimplex::isFixed:
2819                    break;
2820               case isFree:
2821               case superBasic:
2822                    break;
2823               case atUpperBound:
2824                    if (fabs(value - upperValue) > primalTolerance_)
2825                         numberInfeasibilities++;
2826                    break;
2827               case atLowerBound:
2828                    if (fabs(value - lowerValue) > primalTolerance_)
2829                         numberInfeasibilities++;
2830                    break;
2831               }
2832          }
2833          // If dual infeasible then carry on
2834          if (numberInfeasibilities) {
2835               handler_->message(CLP_DUAL_CHECKB, messages_)
2836                         << newBound
2837                         << CoinMessageEol;
2838               int iSequence;
2839               for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
2840                    double lowerValue = lower_[iSequence];
2841                    double upperValue = upper_[iSequence];
2842                    double newLowerValue;
2843                    double newUpperValue;
2844                    Status status = getStatus(iSequence);
2845                    if (status == atUpperBound ||
2846                              status == atLowerBound) {
2847                         double value = solution_[iSequence];
2848                         if (value - lowerValue <= upperValue - value) {
2849                              newLowerValue = CoinMax(lowerValue, value - 0.666667 * newBound);
2850                              newUpperValue = CoinMin(upperValue, newLowerValue + newBound);
2851                         } else {
2852                              newUpperValue = CoinMin(upperValue, value + 0.666667 * newBound);
2853                              newLowerValue = CoinMax(lowerValue, newUpperValue - newBound);
2854                         }
2855                         lower_[iSequence] = newLowerValue;
2856                         upper_[iSequence] = newUpperValue;
2857                         if (newLowerValue > lowerValue) {
2858                              if (newUpperValue < upperValue) {
2859                                   setFakeBound(iSequence, ClpSimplexDual::bothFake);
2860#ifdef CLP_INVESTIGATE
2861                                   abort(); // No idea what should happen here - I have never got here
2862#endif
2863                                   numberFake_++;
2864                              } else {
2865                                   setFakeBound(iSequence, ClpSimplexDual::lowerFake);
2866                                   numberFake_++;
2867                              }
2868                         } else {
2869                              if (newUpperValue < upperValue) {
2870                                   setFakeBound(iSequence, ClpSimplexDual::upperFake);
2871                                   numberFake_++;
2872                              }
2873                         }
2874                         if (status == atUpperBound)
2875                              solution_[iSequence] = newUpperValue;
2876                         else
2877                              solution_[iSequence] = newLowerValue;
2878                         double movement = solution_[iSequence] - value;
2879                         if (movement && outputArray) {
2880                              if (iSequence >= numberColumns_) {
2882                                   changeCost += movement * cost_[iSequence];
2883                              } else {
2884                                   matrix_->add(this, outputArray, iSequence, movement);
2885                                   changeCost += movement * cost_[iSequence];
2886                              }
2887                         }
2888                    }
2889               }
2890               dualBound_ = newBound;
2891          } else {
2892               numberInfeasibilities = -1;
2893          }
2894          return numberInfeasibilities;
2895     } else if (initialize == 1 || initialize == 3) {
2896          int iSequence;
2897          if (initialize == 3) {
2898               for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
2899                    setFakeBound(iSequence, ClpSimplexDual::noFake);
2900               }
2901          }
2902          double testBound = 0.999999 * dualBound_;
2903          for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
2904               Status status = getStatus(iSequence);
2905               if (status == atUpperBound ||
2906                         status == atLowerBound) {
2907                    double lowerValue = lower_[iSequence];
2908                    double upperValue = upper_[iSequence];
2909                    double value = solution_[iSequence];
2910                    if (lowerValue > -largeValue_ || upperValue < largeValue_) {
2911                         if (true || lowerValue - value > -0.5 * dualBound_ ||
2912                                   upperValue - value < 0.5 * dualBound_) {
2913                              if (fabs(lowerValue - value) <= fabs(upperValue - value)) {
2914                                   if (upperValue > lowerValue + testBound) {
2915                                        if (getFakeBound(iSequence) == ClpSimplexDual::noFake)
2916                                             numberFake_++;
2917                                        upper_[iSequence] = lowerValue + dualBound_;
2918                                        setFakeBound(iSequence, ClpSimplexDual::upperFake);
2919                                   }
2920                              } else {
2921                                   if (lowerValue < upperValue - testBound) {
2922                                        if (getFakeBound(iSequence) == ClpSimplexDual::noFake)
2923                                             numberFake_++;
2924                                        lower_[iSequence] = upperValue - dualBound_;
2925                                        setFakeBound(iSequence, ClpSimplexDual::lowerFake);
2926                                   }
2927                              }
2928                         } else {
2929                              if (getFakeBound(iSequence) == ClpSimplexDual::noFake)
2930                                   numberFake_++;
2931                              lower_[iSequence] = -0.5 * dualBound_;
2932                              upper_[iSequence] = 0.5 * dualBound_;
2933                              setFakeBound(iSequence, ClpSimplexDual::bothFake);
2934                              abort();
2935                         }
2936                         if (status == atUpperBound)
2937                              solution_[iSequence] = upper_[iSequence];
2938                         else
2939                              solution_[iSequence] = lower_[iSequence];
2940                    } else {
2941                         // set non basic free variables to fake bounds
2942                         // I don't think we should ever get here
2943                         CoinAssert(!("should not be here"));
2944                         lower_[iSequence] = -0.5 * dualBound_;
2945                         upper_[iSequence] = 0.5 * dualBound_;
2946                         setFakeBound(iSequence, ClpSimplexDual::bothFake);
2947                         numberFake_++;
2948                         setStatus(iSequence, atUpperBound);
2949                         solution_[iSequence] = 0.5 * dualBound_;
2950                    }
2951               } else if (status == basic) {
2952                    // make sure not at fake bound and bounds correct
2953                    setFakeBound(iSequence, ClpSimplexDual::noFake);
2954                    double gap = upper_[iSequence] - lower_[iSequence];
2955                    if (gap > 0.5 * dualBound_ && gap < 2.0 * dualBound_) {
2956                         if (iSequence < numberColumns_) {
2957                              if (columnScale_) {
2958                                   double multiplier = rhsScale_ * inverseColumnScale_[iSequence];
2959                                   // lower
2960                                   double value = columnLower_[iSequence];
2961                                   if (value > -1.0e30) {
2962                                        value *= multiplier;
2963                                   }
2964                                   lower_[iSequence] = value;
2965                                   // upper
2966                                   value = columnUpper_[iSequence];
2967                                   if (value < 1.0e30) {
2968                                        value *= multiplier;
2969                                   }
2970                                   upper_[iSequence] = value;
2971                              } else {
2972                                   lower_[iSequence] = columnLower_[iSequence];;
2973                                   upper_[iSequence] = columnUpper_[iSequence];;
2974                              }
2975                         } else {
2976                              int iRow = iSequence - numberColumns_;
2977                              if (rowScale_) {
2978                                   // lower
2979                                   double multiplier = rhsScale_ * rowScale_[iRow];
2980                                   double value = rowLower_[iRow];
2981                                   if (value > -1.0e30) {
2982                                        value *= multiplier;
2983                                   }
2984                                   lower_[iSequence] = value;
2985                                   // upper
2986                                   value = rowUpper_[iRow];
2987                                   if (value < 1.0e30) {
2988                                        value *= multiplier;
2989                                   }
2990                                   upper_[iSequence] = value;
2991                              } else {
2992                                   lower_[iSequence] = rowLower_[iRow];;
2993                                   upper_[iSequence] = rowUpper_[iRow];;
2994                              }
2995                         }
2996                    }
2997               }
2998          }
2999
3000          return 1;
3001     } else {
3002          // just reset changed ones
3003          if (columnScale_) {
3004               int iSequence;
3005               for (iSequence = 0; iSequence < numberColumns_; iSequence++) {
3006                    FakeBound fakeStatus = getFakeBound(iSequence);
3007                    if (fakeStatus != noFake) {
3008                         if ((static_cast<int> (fakeStatus) & 1) != 0) {
3009                              // lower
3010                              double value = columnLower_[iSequence];
3011                              if (value > -1.0e30) {
3012                                   double multiplier = rhsScale_ * inverseColumnScale_[iSequence];
3013                                   value *= multiplier;
3014                              }
3015                              columnLowerWork_[iSequence] = value;
3016                         }
3017                         if ((static_cast<int> (fakeStatus) & 2) != 0) {
3018                              // upper
3019                              double value = columnUpper_[iSequence];
3020                              if (value < 1.0e30) {
3021                                   double multiplier = rhsScale_ * inverseColumnScale_[iSequence];
3022                                   value *= multiplier;
3023                              }
3024                              columnUpperWork_[iSequence] = value;
3025                         }
3026                    }
3027               }
3028               for (iSequence = 0; iSequence < numberRows_; iSequence++) {
3029                    FakeBound fakeStatus = getFakeBound(iSequence + numberColumns_);
3030                    if (fakeStatus != noFake) {
3031                         if ((static_cast<int> (fakeStatus) & 1) != 0) {
3032                              // lower
3033                              double value = rowLower_[iSequence];
3034                              if (value > -1.0e30) {
3035                                   double multiplier = rhsScale_ * rowScale_[iSequence];
3036                                   value *= multiplier;
3037                              }
3038                              rowLowerWork_[iSequence] = value;
3039                         }
3040                         if ((static_cast<int> (fakeStatus) & 2) != 0) {
3041                              // upper
3042                              double value = rowUpper_[iSequence];
3043                              if (value < 1.0e30) {
3044                                   double multiplier = rhsScale_ * rowScale_[iSequence];
3045                                   value *= multiplier;
3046                              }
3047                              rowUpperWork_[iSequence] = value;
3048                         }
3049                    }
3050               }
3051          } else {
3052               int iSequence;
3053               for (iSequence = 0; iSequence < numberColumns_; iSequence++) {
3054                    FakeBound fakeStatus = getFakeBound(iSequence);
3055                    if ((static_cast<int> (fakeStatus) & 1) != 0) {
3056                         // lower
3057                         columnLowerWork_[iSequence] = columnLower_[iSequence];
3058                    }
3059                    if ((static_cast<int> (fakeStatus) & 2) != 0) {
3060                         // upper
3061                         columnUpperWork_[iSequence] = columnUpper_[iSequence];
3062                    }
3063               }
3064               for (iSequence = 0; iSequence < numberRows_; iSequence++) {
3065                    FakeBound fakeStatus = getFakeBound(iSequence + numberColumns_);
3066                    if ((static_cast<int> (fakeStatus) & 1) != 0) {
3067                         // lower
3068                         rowLowerWork_[iSequence] = rowLower_[iSequence];
3069                    }
3070                    if ((static_cast<int> (fakeStatus) & 2) != 0) {
3071                         // upper
3072                         rowUpperWork_[iSequence] = rowUpper_[iSequence];
3073                    }
3074               }
3075          }
3076          return 0;
3077     }
3078}
3079int
3080ClpSimplexDual::dualColumn0(const CoinIndexedVector * rowArray,
3081                            const CoinIndexedVector * columnArray,
3082                            CoinIndexedVector * spareArray,
3083                            double acceptablePivot,
3084                            double & upperReturn, double &bestReturn, double & badFree)
3085{
3086     // do first pass to get possibles
3087     double * spare = spareArray->denseVector();
3088     int * index = spareArray->getIndices();
3089     const double * work;
3090     int number;
3091     const int * which;
3092     const double * reducedCost;
3093     // We can also see if infeasible or pivoting on free
3094     double tentativeTheta = 1.0e15;
3095     double upperTheta = 1.0e31;
3096     double freePivot = acceptablePivot;
3097     double bestPossible = 0.0;
3098     int numberRemaining = 0;
3099     int i;
3100     badFree = 0.0;
3101     if ((moreSpecialOptions_ & 8) != 0) {
3102          // No free or super basic
3103          double multiplier[] = { -1.0, 1.0};
3104          double dualT = - dualTolerance_;
3105          for (int iSection = 0; iSection < 2; iSection++) {
3106
3108               unsigned char * statusArray;
3109               if (!iSection) {
3110                    work = rowArray->denseVector();
3111                    number = rowArray->getNumElements();
3112                    which = rowArray->getIndices();
3113                    reducedCost = rowReducedCost_;
3114                    addSequence = numberColumns_;
3115                    statusArray = status_ + numberColumns_;
3116               } else {
3117                    work = columnArray->denseVector();
3118                    number = columnArray->getNumElements();
3119                    which = columnArray->getIndices();
3120                    reducedCost = reducedCostWork_;
3121                    addSequence = 0;
3122                    statusArray = status_;
3123               }
3124
3125               for (i = 0; i < number; i++) {
3126                    int iSequence = which[i];
3127                    double alpha;
3128                    double oldValue;
3129                    double value;
3130
3131                    assert (getStatus(iSequence + addSequence) != isFree
3132                            && getStatus(iSequence + addSequence) != superBasic);
3133                    int iStatus = (statusArray[iSequence] & 3) - 1;
3134                    if (iStatus) {
3135                         double mult = multiplier[iStatus-1];
3136                         alpha = work[i] * mult;
3137                         if (alpha > 0.0) {
3138                              oldValue = reducedCost[iSequence] * mult;
3139                              value = oldValue - tentativeTheta * alpha;
3140                              if (value < dualT) {
3141                                   bestPossible = CoinMax(bestPossible, alpha);
3142                                   value = oldValue - upperTheta * alpha;
3143                                   if (value < dualT && alpha >= acceptablePivot) {
3144                                        upperTheta = (oldValue - dualT) / alpha;
3145                                        //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
3146                                   }
3147                                   // add to list
3148                                   spare[numberRemaining] = alpha * mult;
3149                                   index[numberRemaining++] = iSequence + addSequence;
3150                              }
3151                         }
3152                    }
3153               }
3154          }
3155     } else {
3156          // some free or super basic
3157          for (int iSection = 0; iSection < 2; iSection++) {
3158
3160
3161               if (!iSection) {
3162                    work = rowArray->denseVector();
3163                    number = rowArray->getNumElements();
3164                    which = rowArray->getIndices();
3165                    reducedCost = rowReducedCost_;
3166                    addSequence = numberColumns_;
3167               } else {
3168                    work = columnArray->denseVector();
3169                    number = columnArray->getNumElements();
3170                    which = columnArray->getIndices();
3171                    reducedCost = reducedCostWork_;
3172                    addSequence = 0;
3173               }
3174
3175               for (i = 0; i < number; i++) {
3176                    int iSequence = which[i];
3177                    double alpha;
3178                    double oldValue;
3179                    double value;
3180                    bool keep;
3181
3182                    switch(getStatus(iSequence + addSequence)) {
3183
3184                    case basic:
3185                    case ClpSimplex::isFixed:
3186                         break;
3187                    case isFree:
3188                    case superBasic:
3189                         alpha = work[i];
3190                         bestPossible = CoinMax(bestPossible, fabs(alpha));
3191                         oldValue = reducedCost[iSequence];
3192                         // If free has to be very large - should come in via dualRow
3194                         //break;
3195                         if (oldValue > dualTolerance_) {
3196                              keep = true;
3197                         } else if (oldValue < -dualTolerance_) {
3198                              keep = true;
3199                         } else {
3200                              if (fabs(alpha) > CoinMax(10.0 * acceptablePivot, 1.0e-5)) {
3201                                   keep = true;
3202                              } else {
3203                                   keep = false;
3205                              }
3206                         }
3207                         if (keep) {
3208                              // free - choose largest
3209                              if (fabs(alpha) > freePivot) {
3210                                   freePivot = fabs(alpha);
3211                                   sequenceIn_ = iSequence + addSequence;
3212                                   theta_ = oldValue / alpha;
3213                                   alpha_ = alpha;
3214                              }
3215                         }
3216                         break;
3217                    case atUpperBound:
3218                         alpha = work[i];
3219                         oldValue = reducedCost[iSequence];
3220                         value = oldValue - tentativeTheta * alpha;
3221                         //assert (oldValue<=dualTolerance_*1.0001);
3222                         if (value > dualTolerance_) {
3223                              bestPossible = CoinMax(bestPossible, -alpha);
3224                              value = oldValue - upperTheta * alpha;
3225                              if (value > dualTolerance_ && -alpha >= acceptablePivot) {
3226                                   upperTheta = (oldValue - dualTolerance_) / alpha;
3227                                   //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
3228                              }
3229                              // add to list
3230                              spare[numberRemaining] = alpha;
3231                              index[numberRemaining++] = iSequence + addSequence;
3232                         }
3233                         break;
3234                    case atLowerBound:
3235                         alpha = work[i];
3236                         oldValue = reducedCost[iSequence];
3237                         value = oldValue - tentativeTheta * alpha;
3238                         //assert (oldValue>=-dualTolerance_*1.0001);
3239                         if (value < -dualTolerance_) {
3240                              bestPossible = CoinMax(bestPossible, alpha);
3241                              value = oldValue - upperTheta * alpha;
3242                              if (value < -dualTolerance_ && alpha >= acceptablePivot) {
3243                                   upperTheta = (oldValue + dualTolerance_) / alpha;
3244                                   //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
3245                              }
3246                              // add to list
3247                              spare[numberRemaining] = alpha;
3248                              index[numberRemaining++] = iSequence + addSequence;
3249                         }
3250                         break;
3251                    }
3252               }
3253          }
3254     }
3255     upperReturn = upperTheta;
3256     bestReturn = bestPossible;
3257     return numberRemaining;
3258}
3259/*
3260   Row array has row part of pivot row (as duals so sign may be switched)
3261   Column array has column part.
3262   This chooses pivot column.
3263   Spare array will be needed when we start getting clever.
3264   We will check for basic so spare array will never overflow.
3265   If necessary will modify costs
3266*/
3267double
3268ClpSimplexDual::dualColumn(CoinIndexedVector * rowArray,
3269                           CoinIndexedVector * columnArray,
3270                           CoinIndexedVector * spareArray,
3271                           CoinIndexedVector * spareArray2,
3272                           double acceptablePivot,
3273                           CoinBigIndex * /*dubiousWeights*/)
3274{
3275     int numberPossiblySwapped = 0;
3276     int numberRemaining = 0;
3277
3278     double totalThru = 0.0; // for when variables flip
3279     //double saveAcceptable=acceptablePivot;
3280     //acceptablePivot=1.0e-9;
3281
3282     double bestEverPivot = acceptablePivot;
3283     int lastSequence = -1;
3284     double lastPivot = 0.0;
3285     double upperTheta;
3286     double newTolerance = dualTolerance_;
3287     //newTolerance = dualTolerance_+1.0e-6*dblParam_[ClpDualTolerance];
3288     // will we need to increase tolerance
3289     //bool thisIncrease = false;
3290     // If we think we need to modify costs (not if something from broad sweep)
3291     bool modifyCosts = false;
3292     // Increase in objective due to swapping bounds (may be negative)
3293     double increaseInObjective = 0.0;
3294
3295     // use spareArrays to put ones looked at in
3296     // we are going to flip flop between
3297     int iFlip = 0;
3298     // Possible list of pivots
3299     int interesting[2];
3300     // where possible swapped ones are
3301     int swapped[2];
3302     // for zeroing out arrays after
3303     int marker[2][2];
3304     // pivot elements
3305     double * array[2], * spare, * spare2;
3306     // indices
3307     int * indices[2], * index, * index2;
3308     spareArray2->clear();
3309     array[0] = spareArray->denseVector();
3310     indices[0] = spareArray->getIndices();
3311     spare = array[0];
3312     index = indices[0];
3313     array[1] = spareArray2->denseVector();
3314     indices[1] = spareArray2->getIndices();
3315     int i;
3316
3317     // initialize lists
3318     for (i = 0; i < 2; i++) {
3319          interesting[i] = 0;
3320          swapped[i] = numberColumns_;
3321          marker[i][0] = 0;
3322          marker[i][1] = numberColumns_;
3323     }
3324     /*
3325       First we get a list of possible pivots.  We can also see if the
3326       problem looks infeasible or whether we want to pivot in free variable.
3327       This may make objective go backwards but can only happen a finite
3328       number of times and I do want free variables basic.
3329
3330       Then we flip back and forth.  At the start of each iteration
3331       interesting[iFlip] should have possible candidates and swapped[iFlip]
3332       will have pivots if we decide to take a previous pivot.
3333       At end of each iteration interesting[1-iFlip] should have
3334       candidates if we go through this theta and swapped[1-iFlip]
3335       pivots if we don't go through.
3336
3337       At first we increase theta and see what happens.  We start
3338       theta at a reasonable guess.  If in right area then we do bit by bit.
3339
3340      */
3341
3342     // do first pass to get possibles
3343     upperTheta = 1.0e31;
3344     double bestPossible = 0.0;
3345     double badFree = 0.0;
3346     alpha_ = 0.0;
3347     if (spareIntArray_[0] >= 0) {
3348          numberRemaining = dualColumn0(rowArray, columnArray, spareArray,
3349                                        acceptablePivot, upperTheta, bestPossible, badFree);
3350     } else {
3351          // already done
3352          numberRemaining = spareArray->getNumElements();
3353          spareArray->setNumElements(0);
3354          upperTheta = spareDoubleArray_[0];
3355          bestPossible = spareDoubleArray_[1];
3356          if (spareIntArray_[0] == -1) {
3357               theta_ = spareDoubleArray_[2];
3358               alpha_ = spareDoubleArray_[3];
3359               sequenceIn_ = spareIntArray_[1];
3360          } else {
3361#if 0
3362               int n = numberRemaining;
3363               double u = upperTheta;
3364               double b = bestPossible;
3365               upperTheta = 1.0e31;
3366               bestPossible = 0.0;
3367               numberRemaining = dualColumn0(rowArray, columnArray, spareArray,
3368                                             acceptablePivot, upperTheta, bestPossible, badFree);
3369               assert (n == numberRemaining);
3370               assert (fabs(b - bestPossible) < 1.0e-7);
3371               assert (fabs(u - upperTheta) < 1.0e-7);
3372#endif
3373          }
3374     }
3375     // switch off
3376     spareIntArray_[0] = 0;
3377     // We can also see if infeasible or pivoting on free
3378     double tentativeTheta = 1.0e25;
3379     interesting[0] = numberRemaining;
3380     marker[0][0] = numberRemaining;
3381
3382     if (!numberRemaining && sequenceIn_ < 0)
3383          return 0.0; // Looks infeasible
3384
3385     // If sum of bad small pivots too much
3386#define MORE_CAREFUL
3387#ifdef MORE_CAREFUL
3388     bool badSumPivots = false;
3389#endif
3390     if (sequenceIn_ >= 0) {
3391          // free variable - always choose
3392     } else {
3393
3394          theta_ = 1.0e50;
3395          // now flip flop between spare arrays until reasonable theta
3396          tentativeTheta = CoinMax(10.0 * upperTheta, 1.0e-7);
3397
3398          // loops increasing tentative theta until can't go through
3399
3400          while (tentativeTheta < 1.0e22) {
3401               double thruThis = 0.0;
3402
3403               double bestPivot = acceptablePivot;
3404               int bestSequence = -1;
3405
3406               numberPossiblySwapped = numberColumns_;
3407               numberRemaining = 0;
3408
3409               upperTheta = 1.0e50;
3410
3411               spare = array[iFlip];
3412               index = indices[iFlip];
3413               spare2 = array[1-iFlip];
3414               index2 = indices[1-iFlip];
3415
3416               // try 3 different ways
3417               // 1 bias increase by ones with slightly wrong djs
3418               // 2 bias by all
3419               // 3 bias by all - tolerance
3420#define TRYBIAS 3
3421
3422
3423               double increaseInThis = 0.0; //objective increase in this loop
3424
3425               for (i = 0; i < interesting[iFlip]; i++) {
3426                    int iSequence = index[i];
3427                    double alpha = spare[i];
3428                    double oldValue = dj_[iSequence];
3429                    double value = oldValue - tentativeTheta * alpha;
3430
3431                    if (alpha < 0.0) {
3432                         //at upper bound
3433                         if (value > newTolerance) {
3434                              double range = upper_[iSequence] - lower_[iSequence];
3435                              thruThis -= range * alpha;
3436#if TRYBIAS==1
3437                              if (oldValue > 0.0)
3438                                   increaseInThis -= oldValue * range;
3439#elif TRYBIAS==2
3440                              increaseInThis -= oldValue * range;
3441#else
3442                              increaseInThis -= (oldValue + dualTolerance_) * range;
3443#endif
3444                              // goes on swapped list (also means candidates if too many)
3445                              spare2[--numberPossiblySwapped] = alpha;
3446                              index2[numberPossiblySwapped] = iSequence;
3447                              if (fabs(alpha) > bestPivot) {
3448                                   bestPivot = fabs(alpha);
3449                                   bestSequence = numberPossiblySwapped;
3450                              }
3451                         } else {
3452                              value = oldValue - upperTheta * alpha;
3453                              if (value > newTolerance && -alpha >= acceptablePivot)
3454                                   upperTheta = (oldValue - newTolerance) / alpha;
3455                              spare2[numberRemaining] = alpha;
3456                              index2[numberRemaining++] = iSequence;
3457                         }
3458                    } else {
3459                         // at lower bound
3460                         if (value < -newTolerance) {
3461                              double range = upper_[iSequence] - lower_[iSequence];
3462                              thruThis += range * alpha;
3463                              //?? is this correct - and should we look at good ones
3464#if TRYBIAS==1
3465                              if (oldValue < 0.0)
3466                                   increaseInThis += oldValue * range;
3467#elif TRYBIAS==2
3468                              increaseInThis += oldValue * range;
3469#else
3470                              increaseInThis += (oldValue - dualTolerance_) * range;
3471#endif
3472                              // goes on swapped list (also means candidates if too many)
3473                              spare2[--numberPossiblySwapped] = alpha;
3474                              index2[numberPossiblySwapped] = iSequence;
3475                              if (fabs(alpha) > bestPivot) {
3476                                   bestPivot = fabs(alpha);
3477                                   bestSequence = numberPossiblySwapped;
3478                              }
3479                         } else {
3480                              value = oldValue - upperTheta * alpha;
3481                              if (value < -newTolerance && alpha >= acceptablePivot)
3482                                   upperTheta = (oldValue + newTolerance) / alpha;
3483                              spare2[numberRemaining] = alpha;
3484                              index2[numberRemaining++] = iSequence;
3485                         }
3486                    }
3487               }
3488               swapped[1-iFlip] = numberPossiblySwapped;
3489               interesting[1-iFlip] = numberRemaining;
3490               marker[1-iFlip][0] = CoinMax(marker[1-iFlip][0], numberRemaining);
3491               marker[1-iFlip][1] = CoinMin(marker[1-iFlip][1], numberPossiblySwapped);
3492
3493               if (totalThru + thruThis >= fabs(dualOut_) ||
3494                         increaseInObjective + increaseInThis < 0.0) {
3495                    // We should be pivoting in this batch
3496                    // so compress down to this lot
3497                    numberRemaining = 0;
3498                    for (i = numberColumns_ - 1; i >= swapped[1-iFlip]; i--) {
3499                         spare[numberRemaining] = spare2[i];
3500                         index[numberRemaining++] = index2[i];
3501                    }
3502                    interesting[iFlip] = numberRemaining;
3503                    int iTry;
3504#define MAXTRY 100
3505                    // first get ratio with tolerance
3506                    for (iTry = 0; iTry < MAXTRY; iTry++) {
3507
3508                         upperTheta = 1.0e50;
3509                         numberPossiblySwapped = numberColumns_;
3510                         numberRemaining = 0;
3511
3512                         increaseInThis = 0.0; //objective increase in this loop
3513
3514                         thruThis = 0.0;
3515
3516                         spare = array[iFlip];
3517                         index = indices[iFlip];
3518                         spare2 = array[1-iFlip];
3519                         index2 = indices[1-iFlip];
3520                         for (i = 0; i < interesting[iFlip]; i++) {
3521                              int iSequence = index[i];
3522                              double alpha = spare[i];
3523                              double oldValue = dj_[iSequence];
3524                              double value = oldValue - upperTheta * alpha;
3525
3526                              if (alpha < 0.0) {
3527                                   //at upper bound
3528                                   if (value > newTolerance) {
3529                                        if (-alpha >= acceptablePivot) {
3530                                             upperTheta = (oldValue - newTolerance) / alpha;
3531                                        }
3532                                   }
3533                              } else {
3534                                   // at lower bound
3535                                   if (value < -newTolerance) {
3536                                        if (alpha >= acceptablePivot) {
3537                                             upperTheta = (oldValue + newTolerance) / alpha;
3538                                        }
3539                                   }
3540                              }
3541                         }
3542                         bestPivot = acceptablePivot;
3543                         sequenceIn_ = -1;
3544#ifdef DUBIOUS_WEIGHTS
3545                         double bestWeight = COIN_DBL_MAX;
3546#endif
3547                         double largestPivot = acceptablePivot;
3548                         // now choose largest and sum all ones which will go through
3549                         //printf("XX it %d number %d\n",numberIterations_,interesting[iFlip]);
3550                         // Sum of bad small pivots
3551#ifdef MORE_CAREFUL
3552                         double sumBadPivots = 0.0;
3553                         badSumPivots = false;
3554#endif
3555                         // Make sure upperTheta will work (-O2 and above gives problems)
3556                         upperTheta *= 1.0000000001;
3557                         for (i = 0; i < interesting[iFlip]; i++) {
3558                              int iSequence = index[i];
3559                              double alpha = spare[i];
3560                              double value = dj_[iSequence] - upperTheta * alpha;
3561                              double badDj = 0.0;
3562
3563                              bool addToSwapped = false;
3564
3565                              if (alpha < 0.0) {
3566                                   //at upper bound
3567                                   if (value >= 0.0) {
3568                                        addToSwapped = true;
3569#if TRYBIAS==1
3570                                        badDj = -CoinMax(dj_[iSequence], 0.0);
3571#elif TRYBIAS==2
3572                                        badDj = -dj_[iSequence];
3573#else
3574                                        badDj = -dj_[iSequence] - dualTolerance_;
3575#endif
3576                                   }
3577                              } else {
3578                                   // at lower bound
3579                                   if (value <= 0.0) {
3580                                        addToSwapped = true;
3581#if TRYBIAS==1
3582                                        badDj = CoinMin(dj_[iSequence], 0.0);
3583#elif TRYBIAS==2
3584                                        badDj = dj_[iSequence];
3585#else
3586                                        badDj = dj_[iSequence] - dualTolerance_;
3587#endif
3588                                   }
3589                              }
3590                              if (!addToSwapped) {
3591                                   // add to list of remaining
3592                                   spare2[numberRemaining] = alpha;
3593                                   index2[numberRemaining++] = iSequence;
3594                              } else {
3595                                   // add to list of swapped
3596                                   spare2[--numberPossiblySwapped] = alpha;
3597                                   index2[numberPossiblySwapped] = iSequence;
3598                                   // select if largest pivot
3599                                   bool take = false;
3600                                   double absAlpha = fabs(alpha);
3601#ifdef DUBIOUS_WEIGHTS
3602                                   // User could do anything to break ties here
3603                                   double weight;
3604                                   if (dubiousWeights)
3605                                        weight = dubiousWeights[iSequence];
3606                                   else
3607                                        weight = 1.0;
3608                                   weight += randomNumberGenerator_.randomDouble() * 1.0e-2;
3609                                   if (absAlpha > 2.0 * bestPivot) {
3610                                        take = true;
3611                                   } else if (absAlpha > largestPivot) {
3612                                        // could multiply absAlpha and weight
3613                                        if (weight * bestPivot < bestWeight * absAlpha)
3614                                             take = true;
3615                                   }
3616#else
3617                                   if (absAlpha > bestPivot)
3618                                        take = true;
3619#endif
3620#ifdef MORE_CAREFUL
3621                                   if (absAlpha < acceptablePivot && upperTheta < 1.0e20) {
3622                                        if (alpha < 0.0) {
3623                                             //at upper bound
3624                                             if (value > dualTolerance_) {
3625                                                  double gap = upper_[iSequence] - lower_[iSequence];
3626                                                  if (gap < 1.0e20)
3627                                                       sumBadPivots += value * gap;
3628                                                  else
3629                                                       sumBadPivots += 1.0e20;
3630                                                  //printf("bad %d alpha %g dj at upper %g\n",
3631                                                  //     iSequence,alpha,value);
3632                                             }
3633                                        } else {
3634                                             //at lower bound
3635                                             if (value < -dualTolerance_) {
3636                                                  double gap = upper_[iSequence] - lower_[iSequence];
3637                                                  if (gap < 1.0e20)
3638                                                       sumBadPivots -= value * gap;
3639                                                  else
3640                                                       sumBadPivots += 1.0e20;
3641                                                  //printf("bad %d alpha %g dj at lower %g\n",
3642                                                  //     iSequence,alpha,value);
3643                                             }
3644                                        }
3645                                   }
3646#endif
3647#ifdef FORCE_FOLLOW
3648                                   if (iSequence == force_in) {
3649                                        printf("taking %d - alpha %g best %g\n", force_in, absAlpha, largestPivot);
3650                                        take = true;
3651                                   }
3652#endif
3653                                   if (take) {
3654                                        sequenceIn_ = numberPossiblySwapped;
3655                                        bestPivot =  absAlpha;
3656                                        theta_ = dj_[iSequence] / alpha;
3657                                        largestPivot = CoinMax(largestPivot, 0.5 * bestPivot);
3658#ifdef DUBIOUS_WEIGHTS
3659                                        bestWeight = weight;
3660#endif
3661                                        //printf(" taken seq %d alpha %g weight %d\n",
3662                                        //   iSequence,absAlpha,dubiousWeights[iSequence]);
3663                                   } else {
3664                                        //printf(" not taken seq %d alpha %g weight %d\n",
3665                                        //   iSequence,absAlpha,dubiousWeights[iSequence]);
3666                                   }
3667                                   double range = upper_[iSequence] - lower_[iSequence];
3668                                   thruThis += range * fabs(alpha);
3669                                   increaseInThis += badDj * range;
3670                              }
3671                         }
3672                         marker[1-iFlip][0] = CoinMax(marker[1-iFlip][0], numberRemaining);
3673                         marker[1-iFlip][1] = CoinMin(marker[1-iFlip][1], numberPossiblySwapped);
3674#ifdef MORE_CAREFUL
3675                         // If we have done pivots and things look bad set alpha_ 0.0 to force factorization
3676                         if (sumBadPivots > 1.0e4) {
3677                              if (handler_->logLevel() > 1)
3678                                   *handler_ << "maybe forcing re-factorization - sum " << sumBadPivots << " " << factorization_->pivots() << " pivots" << CoinMessageEol;
3679                              if(factorization_->pivots() > 3) {
3680                                   badSumPivots = true;
3681                                   break;
3682                              }
3683                         }
3684#endif
3685                         swapped[1-iFlip] = numberPossiblySwapped;
3686                         interesting[1-iFlip] = numberRemaining;
3687                         // If we stop now this will be increase in objective (I think)
3688                         double increase = (fabs(dualOut_) - totalThru) * theta_;
3689                         increase += increaseInObjective;
3690                         if (theta_ < 0.0)
3691                              thruThis += fabs(dualOut_); // force using this one
3692                         if (increaseInObjective < 0.0 && increase < 0.0 && lastSequence >= 0) {
3693                              // back
3694                              // We may need to be more careful - we could do by
3695                              // switch so we always do fine grained?
3696                              bestPivot = 0.0;
3697                         } else {
3698                              // add in
3699                              totalThru += thruThis;
3700                              increaseInObjective += increaseInThis;
3701                         }
3702                         if (bestPivot < 0.1 * bestEverPivot &&
3703                                   bestEverPivot > 1.0e-6 &&
3704                                   (bestPivot < 1.0e-3 || totalThru * 2.0 > fabs(dualOut_))) {
3705                              // back to previous one
3706                              sequenceIn_ = lastSequence;
3707                              // swap regions
3708                              iFlip = 1 - iFlip;
3709                              break;
3710                         } else if (sequenceIn_ == -1 && upperTheta > largeValue_) {
3711                              if (lastPivot > acceptablePivot) {
3712                                   // back to previous one
3713                                   sequenceIn_ = lastSequence;
3714                                   // swap regions
3715                                   iFlip = 1 - iFlip;
3716                              } else {
3717                                   // can only get here if all pivots too small
3718                              }
3719                              break;
3720                         } else if (totalThru >= fabs(dualOut_)) {
3721                              modifyCosts = true; // fine grain - we can modify costs
3722                              break; // no point trying another loop
3723                         } else {
3724                              lastSequence = sequenceIn_;
3725                              if (bestPivot > bestEverPivot)
3726                                   bestEverPivot = bestPivot;
3727                              iFlip = 1 - iFlip;
3728                              modifyCosts = true; // fine grain - we can modify costs
3729                         }
3730                    }
3731                    if (iTry == MAXTRY)
3732                         iFlip = 1 - iFlip; // flip back
3733                    break;
3734               } else {
3735                    // skip this lot
3736                    if (bestPivot > 1.0e-3 || bestPivot > bestEverPivot) {
3737                         bestEverPivot = bestPivot;
3738                         lastSequence = bestSequence;
3739                    } else {
3740                         // keep old swapped
3741                         CoinMemcpyN(array[iFlip] + swapped[iFlip],
3742                                     numberColumns_ - swapped[iFlip], array[1-iFlip] + swapped[iFlip]);
3743                         CoinMemcpyN(indices[iFlip] + swapped[iFlip],
3744                                     numberColumns_ - swapped[iFlip], indices[1-iFlip] + swapped[iFlip]);
3745                         marker[1-iFlip][1] = CoinMin(marker[1-iFlip][1], swapped[iFlip]);
3746                         swapped[1-iFlip] = swapped[iFlip];
3747                    }
3748                    increaseInObjective += increaseInThis;
3749                    iFlip = 1 - iFlip; // swap regions
3750                    tentativeTheta = 2.0 * upperTheta;
3751                    totalThru += thruThis;
3752               }
3753          }
3754
3755          // can get here without sequenceIn_ set but with lastSequence
3756          if (sequenceIn_ < 0 && lastSequence >= 0) {
3757               // back to previous one
3758               sequenceIn_ = lastSequence;
3759               // swap regions
3760               iFlip = 1 - iFlip;
3761          }
3762
3763#define MINIMUMTHETA 1.0e-18
3764          // Movement should be minimum for anti-degeneracy - unless
3765          // fixed variable out
3766          double minimumTheta;
3767          if (upperOut_ > lowerOut_)
3768               minimumTheta = MINIMUMTHETA;
3769          else
3770               minimumTheta = 0.0;
3771          if (sequenceIn_ >= 0) {
3772               // at this stage sequenceIn_ is just pointer into index array
3773               // flip just so we can use iFlip
3774               iFlip = 1 - iFlip;
3775               spare = array[iFlip];
3776               index = indices[iFlip];
3777               double oldValue;
3778               alpha_ = spare[sequenceIn_];
3779               sequenceIn_ = indices[iFlip][sequenceIn_];
3780               oldValue = dj_[sequenceIn_];
3781               theta_ = CoinMax(oldValue / alpha_, 0.0);
3782               if (theta_ < minimumTheta && fabs(alpha_) < 1.0e5 && 1) {
3783                    // can't pivot to zero
3784#if 0
3785                    if (oldValue - minimumTheta*alpha_ >= -dualTolerance_) {
3786                         theta_ = minimumTheta;
3787                    } else if (oldValue - minimumTheta*alpha_ >= -newTolerance) {
3788                         theta_ = minimumTheta;
3789                         thisIncrease = true;
3790                    } else {
3791                         theta_ = CoinMax((oldValue + newTolerance) / alpha_, 0.0);
3792                         thisIncrease = true;
3793                    }
3794#else
3795                    theta_ = minimumTheta;
3796#endif
3797               }
3798               // may need to adjust costs so all dual feasible AND pivoted is exactly 0
3799               //int costOffset = numberRows_+numberColumns_;
3800               if (modifyCosts && !badSumPivots) {
3801                    int i;
3802                    for (i = numberColumns_ - 1; i >= swapped[iFlip]; i--) {
3803                         int iSequence = index[i];
3804                         double alpha = spare[i];
3805                         double value = dj_[iSequence] - theta_ * alpha;
3806
3807                         // can't be free here
3808
3809                         if (alpha < 0.0) {
3810                              //at upper bound
3811                              if (value > dualTolerance_) {
3812                                   //thisIncrease = true;
3813#if CLP_CAN_HAVE_ZERO_OBJ<2
3814#define MODIFYCOST 2
3815#endif
3816#if MODIFYCOST
3817                                   // modify cost to hit new tolerance
3818                                   double modification = alpha * theta_ - dj_[iSequence]
3819                                                         + newTolerance;
3820                                   if ((specialOptions_&(2048 + 4096 + 16384)) != 0) {
3821                                        if ((specialOptions_ & 16384) != 0) {
3822                                             if (fabs(modification) < 1.0e-8)
3823                                                  modification = 0.0;
3824                                        } else if ((specialOptions_ & 2048) != 0) {
3825                                             if (fabs(modification) < 1.0e-10)
3826                                                  modification = 0.0;
3827                                        } else {
3828                                             if (fabs(modification) < 1.0e-12)
3829                                                  modification = 0.0;
3830                                        }
3831                                   }
3832                                   dj_[iSequence] += modification;
3833                                   cost_[iSequence] +=  modification;
3834                                   if (modification)
3835                                        numberChanged_ ++; // Say changed costs
3836                                   //cost_[iSequence+costOffset] += modification; // save change
3837#endif
3838                              }
3839                         } else {
3840                              // at lower bound
3841                              if (-value > dualTolerance_) {
3842                                   //thisIncrease = true;
3843#if MODIFYCOST
3844                                   // modify cost to hit new tolerance
3845                                   double modification = alpha * theta_ - dj_[iSequence]
3846                                                         - newTolerance;
3847                                   //modification = CoinMax(modification,-dualTolerance_);
3848                                   //assert (fabs(modification)<1.0e-7);
3849                                   if ((specialOptions_&(2048 + 4096)) != 0) {
3850                                        if ((specialOptions_ & 2048) != 0) {
3851                                             if (fabs(modification) < 1.0e-10)
3852                                                  modification = 0.0;
3853                                        } else {
3854                                             if (fabs(modification) < 1.0e-12)
3855                                                  modification = 0.0;
3856                                        }
3857                                   }
3858                                   dj_[iSequence] += modification;
3859                                   cost_[iSequence] +=  modification;
3860                                   if (modification)
3861                                        numberChanged_ ++; // Say changed costs
3862                                   //cost_[iSequence+costOffset] += modification; // save change
3863#endif
3864                              }
3865                         }
3866                    }
3867               }
3868          }
3869     }
3870
3871#ifdef MORE_CAREFUL
3872     // If we have done pivots and things look bad set alpha_ 0.0 to force factorization
3873     if ((badSumPivots ||
3874               fabs(theta_ * badFree) > 10.0 * dualTolerance_) && factorization_->pivots()) {
3875          if (handler_->logLevel() > 1)
3876               *handler_ << "forcing re-factorization" << CoinMessageEol;
3877          sequenceIn_ = -1;
3878     }
3879#endif
3880     if (sequenceIn_ >= 0) {
3881          lowerIn_ = lower_[sequenceIn_];
3882          upperIn_ = upper_[sequenceIn_];
3883          valueIn_ = solution_[sequenceIn_];
3884          dualIn_ = dj_[sequenceIn_];
3885
3886          if (numberTimesOptimal_) {
3887               // can we adjust cost back closer to original
3888               //*** add coding
3889          }
3890#if MODIFYCOST>1
3891          // modify cost to hit zero exactly
3892          // so (dualIn_+modification)==theta_*alpha_
3893          double modification = theta_ * alpha_ - dualIn_;
3894          // But should not move objective too much ??
3895#define DONT_MOVE_OBJECTIVE
3896#ifdef DONT_MOVE_OBJECTIVE
3897          double moveObjective = fabs(modification * solution_[sequenceIn_]);
3898          double smallMove = CoinMax(fabs(objectiveValue_), 1.0e-3);
3899          if (moveObjective > smallMove) {
3900               if (handler_->logLevel() > 1)
3901                    printf("would move objective by %g - original mod %g sol value %g\n", moveObjective,
3902                           modification, solution_[sequenceIn_]);
3903               modification *= smallMove / moveObjective;
3904          }
3905#endif
3907               modification = 0.0;
3908          if ((specialOptions_&(2048 + 4096)) != 0) {
3909               if ((specialOptions_ & 16384) != 0) {
3910                    // in fast dual
3911                    if (fabs(modification) < 1.0e-7)
3912                         modification = 0.0;
3913               } else if ((specialOptions_ & 2048) != 0) {
3914                    if (fabs(modification) < 1.0e-10)
3915                         modification = 0.0;
3916               } else {
3917                    if (fabs(modification) < 1.0e-12)
3918                         modification = 0.0;
3919               }
3920          }
3921          dualIn_ += modification;
3922          dj_[sequenceIn_] = dualIn_;
3923          cost_[sequenceIn_] += modification;
3924          if (modification)
3925               numberChanged_ ++; // Say changed costs
3926          //int costOffset = numberRows_+numberColumns_;
3927          //cost_[sequenceIn_+costOffset] += modification; // save change
3928          //assert (fabs(modification)<1.0e-6);
3929#ifdef CLP_DEBUG
3930          if ((handler_->logLevel() & 32) && fabs(modification) > 1.0e-15)
3931               printf("exact %d new cost %g, change %g\n", sequenceIn_,
3932                      cost_[sequenceIn_], modification);
3933#endif
3934#endif
3935
3936          if (alpha_ < 0.0) {
3937               // as if from upper bound
3938               directionIn_ = -1;
3939               upperIn_ = valueIn_;
3940          } else {
3941               // as if from lower bound
3942               directionIn_ = 1;
3943               lowerIn_ = valueIn_;
3944          }
3945     } else {
3946          // no pivot
3947          bestPossible = 0.0;
3948          alpha_ = 0.0;
3949     }
3950     //if (thisIncrease)
3951     //dualTolerance_+= 1.0e-6*dblParam_[ClpDualTolerance];
3952
3953     // clear arrays
3954
3955     for (i = 0; i < 2; i++) {
3956          CoinZeroN(array[i], marker[i][0]);
3957          CoinZeroN(array[i] + marker[i][1], numberColumns_ - marker[i][1]);
3958     }
3959     return bestPossible;
3960}
3961#ifdef CLP_ALL_ONE_FILE
3962#undef MAXTRY
3963#endif
3964/* Checks if tentative optimal actually means unbounded
3965   Returns -3 if not, 2 if is unbounded */
3966int
3967ClpSimplexDual::checkUnbounded(CoinIndexedVector * ray,
3968                               CoinIndexedVector * spare,
3969                               double changeCost)
3970{
3971     int status = 2; // say unbounded
3972     factorization_->updateColumn(spare, ray);
3973     // get reduced cost
3974     int i;
3975     int number = ray->getNumElements();
3976     int * index = ray->getIndices();
3977     double * array = ray->denseVector();
3978     for (i = 0; i < number; i++) {
3979          int iRow = index[i];
3980          int iPivot = pivotVariable_[iRow];
3981          changeCost -= cost(iPivot) * array[iRow];
3982     }
3983     double way;
3984     if (changeCost > 0.0) {
3985          //try going down
3986          way = 1.0;
3987     } else if (changeCost < 0.0) {
3988          //try going up
3989          way = -1.0;
3990     } else {
3991#ifdef CLP_DEBUG
3992          printf("can't decide on up or down\n");
3993#endif
3994          way = 0.0;
3995          status = -3;
3996     }
3997     double movement = 1.0e10 * way; // some largish number
3998     double zeroTolerance = 1.0e-14 * dualBound_;
3999     for (i = 0; i < number; i++) {
4000          int iRow = index[i];
4001          int iPivot = pivotVariable_[iRow];
4002          double arrayValue = array[iRow];
4003          if (fabs(arrayValue) < zeroTolerance)
4004               arrayValue = 0.0;
4005          double newValue = solution(iPivot) + movement * arrayValue;
4006          if (newValue > upper(iPivot) + primalTolerance_ ||
4007                    newValue < lower(iPivot) - primalTolerance_)
4008               status = -3; // not unbounded
4009     }
4010     if (status == 2) {
4011          // create ray
4012          delete [] ray_;
4013          ray_ = new double [numberColumns_];
4014          CoinZeroN(ray_, numberColumns_);
4015          for (i = 0; i < number; i++) {
4016               int iRow = index[i];
4017               int iPivot = pivotVariable_[iRow];
4018               double arrayValue = array[iRow];
4019               if (iPivot < numberColumns_ && fabs(arrayValue) >= zeroTolerance)
4020                    ray_[iPivot] = way * array[iRow];
4021          }
4022     }
4023     ray->clear();
4024     return status;
4025}
4026//static int count_alpha=0;
4027/* Checks if finished.  Updates status */
4028void
4029ClpSimplexDual::statusOfProblemInDual(int & lastCleaned, int type,
4030                                      double * givenDuals, ClpDataSave & saveData,
4031                                      int ifValuesPass)
4032{
4033#ifdef CLP_INVESTIGATE_SERIAL
4034     if (z_thinks > 0 && z_thinks < 2)
4035          z_thinks += 2;
4036#endif
4037     bool arraysNotCreated = (type==0);
4038     // If lots of iterations then adjust costs if large ones
4039     if (numberIterations_ > 4 * (numberRows_ + numberColumns_) && objectiveScale_ == 1.0) {
4040          double largest = 0.0;
4041          for (int i = 0; i < numberRows_; i++) {
4042               int iColumn = pivotVariable_[i];
4043               largest = CoinMax(largest, fabs(cost_[iColumn]));
4044          }
4045          if (largest > 1.0e6) {
4046               objectiveScale_ = 1.0e6 / largest;
4047               for (int i = 0; i < numberRows_ + numberColumns_; i++)
4048                    cost_[i] *= objectiveScale_;
4049          }
4050     }
4051     int numberPivots = factorization_->pivots();
4052     double realDualInfeasibilities = 0.0;
4053     if (type == 2) {
4054          if (alphaAccuracy_ != -1.0)
4055               alphaAccuracy_ = -2.0;
4056          // trouble - restore solution
4057          CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4058          CoinMemcpyN(savedSolution_ + numberColumns_ ,
4059                      numberRows_, rowActivityWork_);
4060          CoinMemcpyN(savedSolution_ ,
4061                      numberColumns_, columnActivityWork_);
4062          // restore extra stuff
4063          int dummy;
4064          matrix_->generalExpanded(this, 6, dummy);
4065          forceFactorization_ = 1; // a bit drastic but ..
4066          changeMade_++; // say something changed
4067          // get correct bounds on all variables
4068          resetFakeBounds(0);
4069     }
4070     int tentativeStatus = problemStatus_;
4071     double changeCost;
4072     bool unflagVariables = true;
4073     bool weightsSaved = false;
4074     bool weightsSaved2 = numberIterations_ && !numberPrimalInfeasibilities_;
4075     int dontFactorizePivots = dontFactorizePivots_;
4076     if (type == 3) {
4077          type = 1;
4078          dontFactorizePivots = 1;
4079     }
4080     if (alphaAccuracy_ < 0.0 || !numberPivots || alphaAccuracy_ > 1.0e4 || numberPivots > 20) {
4081          if (problemStatus_ > -3 || numberPivots > dontFactorizePivots) {
4082               // factorize
4083               // later on we will need to recover from singularities
4084               // also we could skip if first time
4085               // save dual weights
4086               dualRowPivot_->saveWeights(this, 1);
4087               weightsSaved = true;
4088               if (type) {
4089                    // is factorization okay?
4090                    if (internalFactorize(1)) {
4091                         // no - restore previous basis
4092                         unflagVariables = false;
4093                         assert (type == 1);
4094                         changeMade_++; // say something changed
4095                         // Keep any flagged variables
4096                         int i;
4097                         for (i = 0; i < numberRows_ + numberColumns_; i++) {
4098                              if (flagged(i))
4099                                   saveStatus_[i] |= 64; //say flagged
4100                         }
4101                         CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4102                         CoinMemcpyN(savedSolution_ + numberColumns_ ,
4103                                     numberRows_, rowActivityWork_);
4104                         CoinMemcpyN(savedSolution_ ,
4105                                     numberColumns_, columnActivityWork_);
4106                         // restore extra stuff
4107                         int dummy;
4108                         matrix_->generalExpanded(this, 6, dummy);
4109                         // get correct bounds on all variables
4110                         resetFakeBounds(1);
4111                         // need to reject something
4112                         char x = isColumn(sequenceOut_) ? 'C' : 'R';
4113                         handler_->message(CLP_SIMPLEX_FLAG, messages_)
4114                                   << x << sequenceWithin(sequenceOut_)
4115                                   << CoinMessageEol;
4116#ifdef COIN_DEVELOP
4117                         printf("flag d\n");
4118#endif
4119                         setFlagged(sequenceOut_);
4121
4122                         // Go to safe
4123                         factorization_->pivotTolerance(0.99);
4124                         forceFactorization_ = 1; // a bit drastic but ..
4125                         type = 2;
4126                         //assert (internalFactorize(1)==0);
4127                         if (internalFactorize(1)) {
4128                              CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4129                              CoinMemcpyN(savedSolution_ + numberColumns_ ,
4130                                          numberRows_, rowActivityWork_);
4131                              CoinMemcpyN(savedSolution_ ,
4132                                          numberColumns_, columnActivityWork_);
4133                              // restore extra stuff
4134                              int dummy;
4135                              matrix_->generalExpanded(this, 6, dummy);
4136                              // debug
4137                              int returnCode = internalFactorize(1);
4138                              while (returnCode) {
4139                                   // ouch
4140                                   // switch off dense
4141                                   int saveDense = factorization_->denseThreshold();
4142                                   factorization_->setDenseThreshold(0);
4143                                   // Go to safe
4144                                   factorization_->pivotTolerance(0.99);
4145                                   // make sure will do safe factorization
4146                                   pivotVariable_[0] = -1;
4147                                   returnCode = internalFactorize(2);
4148                                   factorization_->setDenseThreshold(saveDense);
4149                              }
4150                              // get correct bounds on all variables
4151                              resetFakeBounds(1);
4152                         }
4153                    }
4154               }
4155               if (problemStatus_ != -4 || numberPivots > 10)
4156                    problemStatus_ = -3;
4157          }
4158     } else {
4159          //printf("testing with accuracy of %g and status of %d\n",alphaAccuracy_,problemStatus_);
4160          //count_alpha++;
4161          //if ((count_alpha%5000)==0)
4162          //printf("count alpha %d\n",count_alpha);
4163     }
4164     // at this stage status is -3 or -4 if looks infeasible
4165     // get primal and dual solutions
4166#if 0
4167     {
4168          int numberTotal = numberRows_ + numberColumns_;
4169          double * saveSol = CoinCopyOfArray(solution_, numberTotal);
4170          double * saveDj = CoinCopyOfArray(dj_, numberTotal);
4171          double tolerance = type ? 1.0e-4 : 1.0e-8;
4172          // always if values pass
4173          double saveObj = objectiveValue_;
4174          double sumPrimal = sumPrimalInfeasibilities_;
4175          int numberPrimal = numberPrimalInfeasibilities_;
4176          double sumDual = sumDualInfeasibilities_;
4177          int numberDual = numberDualInfeasibilities_;
4178          gutsOfSolution(givenDuals, NULL);
4179          int j;
4180          double largestPrimal = tolerance;
4181          int iPrimal = -1;
4182          for (j = 0; j < numberTotal; j++) {
4183               double difference = solution_[j] - saveSol[j];
4184               if (fabs(difference) > largestPrimal) {
4185                    iPrimal = j;
4186                    largestPrimal = fabs(difference);
4187               }
4188          }
4189          double largestDual = tolerance;
4190          int iDual = -1;
4191          for (j = 0; j < numberTotal; j++) {
4192               double difference = dj_[j] - saveDj[j];
4193               if (fabs(difference) > largestDual && upper_[j] > lower_[j]) {
4194                    iDual = j;
4195                    largestDual = fabs(difference);
4196               }
4197          }
4198          if (!type) {
4199               if (fabs(saveObj - objectiveValue_) > 1.0e-5 ||
4200                         numberPrimal != numberPrimalInfeasibilities_ || numberPrimal != 1 ||
4201                         fabs(sumPrimal - sumPrimalInfeasibilities_) > 1.0e-5 || iPrimal >= 0 ||
4202                         numberDual != numberDualInfeasibilities_ || numberDual != 0 ||
4203                         fabs(sumDual - sumDualInfeasibilities_) > 1.0e-5 || iDual >= 0)
4204                    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",
4205                           type, numberIterations_, numberPivots,
4206                           numberPrimal, numberPrimalInfeasibilities_, sumPrimal, sumPrimalInfeasibilities_,
4207                           largestPrimal, iPrimal,
4208                           numberDual, numberDualInfeasibilities_, sumDual, sumDualInfeasibilities_,
4209                           largestDual, iDual,
4210                           saveObj, objectiveValue_);
4211          } else {
4212               if (fabs(saveObj - objectiveValue_) > 1.0e-5 ||
4213                         numberPrimalInfeasibilities_ || iPrimal >= 0 ||
4214                         numberDualInfeasibilities_ || iDual >= 0)
4215                    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",
4216                           type, numberIterations_, numberPivots,
4217                           numberPrimal, numberPrimalInfeasibilities_, sumPrimal, sumPrimalInfeasibilities_,
4218                           largestPrimal, iPrimal,
4219                           numberDual, numberDualInfeasibilities_, sumDual, sumDualInfeasibilities_,
4220                           largestDual, iDual,
4221                           saveObj, objectiveValue_);
4222          }
4223          delete [] saveSol;
4224          delete [] saveDj;
4225     }
4226#else
4227     if (type || ifValuesPass)
4228          gutsOfSolution(givenDuals, NULL);
4229#endif
4230     // If bad accuracy treat as singular
4231     if ((largestPrimalError_ > 1.0e15 || largestDualError_ > 1.0e15) && numberIterations_) {
4232          // restore previous basis
4233          unflagVariables = false;
4234          changeMade_++; // say something changed
4235          // Keep any flagged variables
4236          int i;
4237          for (i = 0; i < numberRows_ + numberColumns_; i++) {
4238               if (flagged(i))
4239                    saveStatus_[i] |= 64; //say flagged
4240          }
4241          CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4242          CoinMemcpyN(savedSolution_ + numberColumns_ ,
4243                      numberRows_, rowActivityWork_);
4244          CoinMemcpyN(savedSolution_ ,
4245                      numberColumns_, columnActivityWork_);
4246          // restore extra stuff
4247          int dummy;
4248          matrix_->generalExpanded(this, 6, dummy);
4249          // get correct bounds on all variables
4250          resetFakeBounds(1);
4251          // need to reject something
4252          char x = isColumn(sequenceOut_) ? 'C' : 'R';
4253          handler_->message(CLP_SIMPLEX_FLAG, messages_)
4254                    << x << sequenceWithin(sequenceOut_)
4255                    << CoinMessageEol;
4256#ifdef COIN_DEVELOP
4257          printf("flag e\n");
4258#endif
4259          setFlagged(sequenceOut_);
4261
4262          // Go to safer
4263          double newTolerance = CoinMin(1.1 * factorization_->pivotTolerance(), 0.99);
4264          factorization_->pivotTolerance(newTolerance);
4265          forceFactorization_ = 1; // a bit drastic but ..
4266          if (alphaAccuracy_ != -1.0)
4267               alphaAccuracy_ = -2.0;
4268          type = 2;
4269          //assert (internalFactorize(1)==0);
4270          if (internalFactorize(1)) {
4271               CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4272               CoinMemcpyN(savedSolution_ + numberColumns_ ,
4273                           numberRows_, rowActivityWork_);
4274               CoinMemcpyN(savedSolution_ ,
4275                           numberColumns_, columnActivityWork_);
4276               // restore extra stuff
4277               int dummy;
4278               matrix_->generalExpanded(this, 6, dummy);
4279               // debug
4280               int returnCode = internalFactorize(1);
4281               while (returnCode) {
4282                    // ouch
4283                    // switch off dense
4284                    int saveDense = factorization_->denseThreshold();
4285                    factorization_->setDenseThreshold(0);
4286                    // Go to safe
4287                    factorization_->pivotTolerance(0.99);
4288                    // make sure will do safe factorization
4289                    pivotVariable_[0] = -1;
4290                    returnCode = internalFactorize(2);
4291                    factorization_->setDenseThreshold(saveDense);
4292               }
4293               // get correct bounds on all variables
4294               resetFakeBounds(1);
4295          }
4296          // get primal and dual solutions
4297          gutsOfSolution(givenDuals, NULL);
4298     } else if (goodAccuracy()) {
4299          // Can reduce tolerance
4300          double newTolerance = CoinMax(0.99 * factorization_->pivotTolerance(), saveData.pivotTolerance_);
4301          factorization_->pivotTolerance(newTolerance);
4302     }
4303     bestObjectiveValue_ = CoinMax(bestObjectiveValue_,
4304                                   objectiveValue_ - bestPossibleImprovement_);
4305     bool reallyBadProblems = false;
4306     // Double check infeasibility if no action
4307     if (progress_.lastIterationNumber(0) == numberIterations_) {
4308          if (dualRowPivot_->looksOptimal()) {
4309               numberPrimalInfeasibilities_ = 0;
4310               sumPrimalInfeasibilities_ = 0.0;
4311          }
4312#if 1
4313     } else {
4314          double thisObj = objectiveValue_ - bestPossibleImprovement_;
4315#ifdef CLP_INVESTIGATE
4316          assert (bestPossibleImprovement_ > -1000.0 && objectiveValue_ > -1.0e100);
4317          if (bestPossibleImprovement_)
4318               printf("obj %g add in %g -> %g\n", objectiveValue_, bestPossibleImprovement_,
4319                      thisObj);
4320#endif
4321          double lastObj = progress_.lastObjective(0);
4322#ifndef NDEBUG
4323#ifdef COIN_DEVELOP
4324          resetFakeBounds(-1);
4325#endif
4326#endif
4327#ifdef CLP_REPORT_PROGRESS
4328          ixxxxxx++;
4329          if (ixxxxxx >= ixxyyyy - 4 && ixxxxxx <= ixxyyyy) {
4330               char temp[20];
4331               sprintf(temp, "sol%d.out", ixxxxxx);
4332               printf("sol%d.out\n", ixxxxxx);
4333               FILE * fp = fopen(temp, "w");
4334               int nTotal = numberRows_ + numberColumns_;
4335               for (int i = 0; i < nTotal; i++)
4336                    fprintf(fp, "%d %d %g %g %g %g %g\n",
4337                            i, status_[i], lower_[i], solution_[i], upper_[i], cost_[i], dj_[i]);
4338               fclose(fp);
4339          }
4340#endif
4341          if(!ifValuesPass && firstFree_ < 0) {
4342               double testTol = 5.0e-3;
4343               if (progress_.timesFlagged() > 10) {
4344                    testTol *= pow(2.0, progress_.timesFlagged() - 8);
4345               } else if (progress_.timesFlagged() > 5) {
4346                    testTol *= 5.0;
4347               }
4348               if (lastObj > thisObj +
4349                         testTol*(fabs(thisObj) + fabs(lastObj)) + testTol) {
4350                    int maxFactor = factorization_->maximumPivots();
4351                    if ((specialOptions_ & 1048576) == 0) {
4352                         if (progress_.timesFlagged() > 10)
4354                         if (maxFactor > 10 - 9) {
4355#ifdef COIN_DEVELOP
4356                              printf("lastobj %g thisobj %g\n", lastObj, thisObj);
4357#endif
4358                              //if (forceFactorization_<0)
4359                              //forceFactorization_= maxFactor;
4360                              //forceFactorization_ = CoinMax(1,(forceFactorization_>>1));
4361                              if ((progressFlag_ & 4) == 0 && lastObj < thisObj + 1.0e4 &&
4362                                        largestPrimalError_ < 1.0e2) {
4363                                   // Just save costs
4364                                   // save extra copy of cost_
4365                                   int nTotal = numberRows_ + numberColumns_;
4366                                   double * temp = new double [2*nTotal];
4367                                   memcpy(temp, cost_, nTotal * sizeof(double));
4368                                   memcpy(temp + nTotal, cost_, nTotal * sizeof(double));
4369                                   delete [] cost_;
4370                                   cost_ = temp;
4371                                   objectiveWork_ = cost_;
4372                                   rowObjectiveWork_ = cost_ + numberColumns_;
4373                                   progressFlag_ |= 4;
4374                              } else {
4375                                   forceFactorization_ = 1;
4376#ifdef COIN_DEVELOP
4377                                   printf("Reducing factorization frequency - bad backwards\n");
4378#endif
4379#if 1
4380                                   unflagVariables = false;
4381                                   changeMade_++; // say something changed
4382                                   int nTotal = numberRows_ + numberColumns_;
4383                                   CoinMemcpyN(saveStatus_, nTotal, status_);
4384                                   CoinMemcpyN(savedSolution_ + numberColumns_ ,
4385                                               numberRows_, rowActivityWork_);
4386                                   CoinMemcpyN(savedSolution_ ,
4387                                               numberColumns_, columnActivityWork_);
4388                                   if ((progressFlag_ & 4) == 0) {
4389                                        // save extra copy of cost_
4390                                        double * temp = new double [2*nTotal];
4391                                        memcpy(temp, cost_, nTotal * sizeof(double));
4392                                        memcpy(temp + nTotal, cost_, nTotal * sizeof(double));
4393                                        delete [] cost_;
4394                                        cost_ = temp;
4395                                        objectiveWork_ = cost_;
4396                                        rowObjectiveWork_ = cost_ + numberColumns_;
4397                                        progressFlag_ |= 4;
4398                                   } else {
4399                                        memcpy(cost_, cost_ + nTotal, nTotal * sizeof(double));
4400                                   }
4401                                   // restore extra stuff
4402                                   int dummy;
4403                                   matrix_->generalExpanded(this, 6, dummy);
4404                                   double pivotTolerance = factorization_->pivotTolerance();
4405                                   if(pivotTolerance < 0.2)
4406                                        factorization_->pivotTolerance(0.2);
4407                                   else if(progress_.timesFlagged() > 2)
4408                                        factorization_->pivotTolerance(CoinMin(pivotTolerance * 1.1, 0.99));
4409                                   if (alphaAccuracy_ != -1.0)
4410                                        alphaAccuracy_ = -2.0;
4411                                   if (internalFactorize(1)) {
4412                                        CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4413                                        CoinMemcpyN(savedSolution_ + numberColumns_ ,
4414                                                    numberRows_, rowActivityWork_);
4415                                        CoinMemcpyN(savedSolution_ ,
4416                                                    numberColumns_, columnActivityWork_);
4417                                        // restore extra stuff
4418                                        int dummy;
4419                                        matrix_->generalExpanded(this, 6, dummy);
4420                                        // debug
4421                                        int returnCode = internalFactorize(1);
4422                                        while (returnCode) {
4423                                             // ouch
4424                                             // switch off dense
4425                                             int saveDense = factorization_->denseThreshold();
4426                                             factorization_->setDenseThreshold(0);
4427                                             // Go to safe
4428                                             factorization_->pivotTolerance(0.99);
4429                                             // make sure will do safe factorization
4430                                             pivotVariable_[0] = -1;
4431                                             returnCode = internalFactorize(2);
4432                                             factorization_->setDenseThreshold(saveDense);
4433                                        }
4434                                   }
4435                                   resetFakeBounds(0);
4436                                   type = 2; // so will restore weights
4437                                   // get primal and dual solutions
4438                                   gutsOfSolution(givenDuals, NULL);
4439                                   if (numberPivots < 2) {
4440                                        // need to reject something
4441                                        char x = isColumn(sequenceOut_) ? 'C' : 'R';
4442                                        handler_->message(CLP_SIMPLEX_FLAG, messages_)
4443                                                  << x << sequenceWithin(sequenceOut_)
4444                                                  << CoinMessageEol;
4445#ifdef COIN_DEVELOP
4446                                        printf("flag d\n");
4447#endif
4448                                        setFlagged(sequenceOut_);
4450                                        progress_.incrementTimesFlagged();
4451                                   }
4452                                   if (numberPivots < 10)
4453                                        reallyBadProblems = true;
4454#ifdef COIN_DEVELOP
4455                                   printf("obj now %g\n", objectiveValue_);
4456#endif
4457                                   progress_.modifyObjective(objectiveValue_
4458                                                             - bestPossibleImprovement_);
4459#endif
4460                              }
4461                         }
4462                    } else {
4463                         // in fast dual give up
4464#ifdef COIN_DEVELOP
4465                         printf("In fast dual?\n");
4466#endif
4467                         problemStatus_ = 3;
4468                    }
4469               } else if (lastObj < thisObj - 1.0e-5 * CoinMax(fabs(thisObj), fabs(lastObj)) - 1.0e-3) {
4470                    numberTimesOptimal_ = 0;
4471               }
4472          }
4473#endif
4474     }
4475     // Up tolerance if looks a bit odd
4476     if (numberIterations_ > CoinMax(1000, numberRows_ >> 4) && (specialOptions_ & 64) != 0) {
4477          if (sumPrimalInfeasibilities_ && sumPrimalInfeasibilities_ < 1.0e5) {
4478               int backIteration = progress_.lastIterationNumber(CLP_PROGRESS - 1);
4479               if (backIteration > 0 && numberIterations_ - backIteration < 9 * CLP_PROGRESS) {
4480                    if (factorization_->pivotTolerance() < 0.9) {
4481                         // up tolerance
4482                         factorization_->pivotTolerance(CoinMin(factorization_->pivotTolerance() * 1.05 + 0.02, 0.91));
4483                         //printf("tol now %g\n",factorization_->pivotTolerance());
4484                         progress_.clearIterationNumbers();
4485                    }
4486               }
4487          }
4488     }
4489     // Check if looping
4490     int loop;
4491     if (!givenDuals && type != 2)
4492          loop = progress_.looping();
4493     else
4494          loop = -1;
4495     if (progress_.reallyBadTimes() > 10) {
4496          problemStatus_ = 10; // instead - try other algorithm
4497#if COIN_DEVELOP>2
4498          printf("returning at %d\n", __LINE__);
4499#endif
4500     }
4501     int situationChanged = 0;
4502     if (loop >= 0) {
4503          problemStatus_ = loop; //exit if in loop
4504          if (!problemStatus_) {
4505               // declaring victory
4506               numberPrimalInfeasibilities_ = 0;
4507               sumPrimalInfeasibilities_ = 0.0;
4508          } else {
4509               problemStatus_ = 10; // instead - try other algorithm
4510#if COIN_DEVELOP>2
4511               printf("returning at %d\n", __LINE__);
4512#endif
4513          }
4514          return;
4515     } else if (loop < -1) {
4516          // something may have changed
4517          gutsOfSolution(NULL, NULL);
4518          situationChanged = 1;
4519     }
4520     // really for free variables in
4521     if((progressFlag_ & 2) != 0) {
4522          situationChanged = 2;
4523     }
4524     progressFlag_ &= (~3); //reset progress flag
4525     if ((progressFlag_ & 4) != 0) {
4526          // save copy of cost_
4527          int nTotal = numberRows_ + numberColumns_;
4528          memcpy(cost_ + nTotal, cost_, nTotal * sizeof(double));
4529     }
4530     /*if (!numberIterations_&&sumDualInfeasibilities_)
4531       printf("OBJ %g sumPinf %g sumDinf %g\n",
4532        objectiveValue(),sumPrimalInfeasibilities_,
4533        sumDualInfeasibilities_);*/
4534     // mark as having gone optimal if looks like it
4535     if (!numberPrimalInfeasibilities_&&
4536         !numberDualInfeasibilities_)
4537       progressFlag_ |= 8;
4538     if (handler_->detail(CLP_SIMPLEX_STATUS, messages_) < 100) {
4539          handler_->message(CLP_SIMPLEX_STATUS, messages_)
4540                    << numberIterations_ << objectiveValue();
4541          handler_->printing(sumPrimalInfeasibilities_ > 0.0)
4542                    << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
4543          handler_->printing(sumDualInfeasibilities_ > 0.0)
4544                    << sumDualInfeasibilities_ << numberDualInfeasibilities_;
4545          handler_->printing(numberDualInfeasibilitiesWithoutFree_
4546                             < numberDualInfeasibilities_)
4547                    << numberDualInfeasibilitiesWithoutFree_;
4548          handler_->message() << CoinMessageEol;
4549     }
4550#if 0
4551     printf("IT %d %g %g(%d) %g(%d)\n",
4552            numberIterations_, objectiveValue(),
4553            sumPrimalInfeasibilities_, numberPrimalInfeasibilities_,
4554            sumDualInfeasibilities_, numberDualInfeasibilities_);
4555#endif
4556     double approximateObjective = objectiveValue_;
4557#ifdef CLP_REPORT_PROGRESS
4558     if (ixxxxxx >= ixxyyyy - 4 && ixxxxxx <= ixxyyyy) {
4559          char temp[20];
4560          sprintf(temp, "x_sol%d.out", ixxxxxx);
4561          FILE * fp = fopen(temp, "w");
4562          int nTotal = numberRows_ + numberColumns_;
4563          for (int i = 0; i < nTotal; i++)
4564               fprintf(fp, "%d %d %g %g %g %g %g\n",
4565                       i, status_[i], lower_[i], solution_[i], upper_[i], cost_[i], dj_[i]);
4566          fclose(fp);
4567          if (ixxxxxx == ixxyyyy)
4568               exit(6);
4569     }
4570#endif
4571     realDualInfeasibilities = sumDualInfeasibilities_;
4572     double saveTolerance = dualTolerance_;
4573     // If we need to carry on cleaning variables
4574     if (!numberPrimalInfeasibilities_ && (specialOptions_ & 1024) != 0 && CLEAN_FIXED) {
4575          for (int iRow = 0; iRow < numberRows_; iRow++) {
4576               int iPivot = pivotVariable_[iRow];
4577               if (!flagged(iPivot) && pivoted(iPivot)) {
4578                    // carry on
4579                    numberPrimalInfeasibilities_ = -1;
4580                    sumOfRelaxedPrimalInfeasibilities_ = 1.0;
4581                    sumPrimalInfeasibilities_ = 1.0;
4582                    break;
4583               }
4584          }
4585     }
4586     /* If we are primal feasible and any dual infeasibilities are on
4587        free variables then it is better to go to primal */
4588     if (!numberPrimalInfeasibilities_ && !numberDualInfeasibilitiesWithoutFree_ &&
4589               numberDualInfeasibilities_)
4590          problemStatus_ = 10;
4591     // dual bound coming in
4592     //double saveDualBound = dualBound_;
4593     bool needCleanFake = false;
4594     while (problemStatus_ <= -3) {
4595          int cleanDuals = 0;
4596          if (situationChanged != 0)
4597               cleanDuals = 1;
4598          int numberChangedBounds = 0;
4599          int doOriginalTolerance = 0;
4600          if ( lastCleaned == numberIterations_)
4601               doOriginalTolerance = 1;
4602          // check optimal
4603          // give code benefit of doubt
4604          if (sumOfRelaxedDualInfeasibilities_ == 0.0 &&
4605                    sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
4606               // say optimal (with these bounds etc)
4607               numberDualInfeasibilities_ = 0;
4608               sumDualInfeasibilities_ = 0.0;
4609               numberPrimalInfeasibilities_ = 0;
4610               sumPrimalInfeasibilities_ = 0.0;
4611          }
4612          //if (dualFeasible()||problemStatus_==-4||(primalFeasible()&&!numberDualInfeasibilitiesWithoutFree_)) {
4613          if (dualFeasible() || problemStatus_ == -4) {
4614               progress_.modifyObjective(objectiveValue_
4615                                         - bestPossibleImprovement_);
4616#ifdef COIN_DEVELOP
4617               if (sumDualInfeasibilities_ || bestPossibleImprovement_)
4618                    printf("improve %g dualinf %g -> %g\n",
4619                           bestPossibleImprovement_, sumDualInfeasibilities_,
4620                           sumDualInfeasibilities_ * dualBound_);
4621#endif
4622               // see if cutoff reached
4623               double limit = 0.0;
4624               getDblParam(ClpDualObjectiveLimit, limit);
4625#if 0
4626               if(fabs(limit) < 1.0e30 && objectiveValue()*optimizationDirection_ >
4627                         limit + 1.0e-7 + 1.0e-8 * fabs(limit) && !numberAtFakeBound()) {
4628                    //looks infeasible on objective
4629                    if (perturbation_ == 101) {
4630                         cleanDuals = 1;
4631                         // Save costs
4632                         int numberTotal = numberRows_ + numberColumns_;
4633                         double * saveCost = CoinCopyOfArray(cost_, numberTotal);
4634                         // make sure fake bounds are back
4635                         changeBounds(1, NULL, changeCost);
4636                         createRim4(false);
4637                         // make sure duals are current
4638                         computeDuals(givenDuals);
4639                         checkDualSolution();
4640                         if(objectiveValue()*optimizationDirection_ >
4641                                   limit + 1.0e-7 + 1.0e-8 * fabs(limit) && !numberDualInfeasibilities_) {
4642                              perturbation_ = 102; // stop any perturbations
4643                              printf("cutoff test succeeded\n");
4644                         } else {
4645                              printf("cutoff test failed\n");
4646                              // put back
4647                              memcpy(cost_, saveCost, numberTotal * sizeof(double));
4648                              // make sure duals are current
4649                              computeDuals(givenDuals);
4650                              checkDualSolution();
4651                              progress_.modifyObjective(-COIN_DBL_MAX);
4652                              problemStatus_ = -1;
4653                         }
4654                         delete [] saveCost;
4655                    }
4656               }
4657#endif
4658               if (primalFeasible() && !givenDuals) {
4659                    // may be optimal - or may be bounds are wrong
4660                    handler_->message(CLP_DUAL_BOUNDS, messages_)
4661                              << dualBound_
4662                              << CoinMessageEol;
4663                    // save solution in case unbounded
4664                    double * saveColumnSolution = NULL;
4665                    double * saveRowSolution = NULL;
4666                    bool inCbc = (specialOptions_ & (0x01000000 | 16384)) != 0;
4667                    if (!inCbc) {
4668                         saveColumnSolution = CoinCopyOfArray(columnActivityWork_, numberColumns_);
4669                         saveRowSolution = CoinCopyOfArray(rowActivityWork_, numberRows_);
4670                    }
4671                    numberChangedBounds = changeBounds(0, rowArray_[3], changeCost);
4672                    if (numberChangedBounds <= 0 && !numberDualInfeasibilities_) {
4673                         //looks optimal - do we need to reset tolerance
4674                         if (perturbation_ == 101) {
4675                              perturbation_ = 102; // stop any perturbations
4676                              cleanDuals = 1;
4677                              // make sure fake bounds are back
4678                              //computeObjectiveValue();
4679                              changeBounds(1, NULL, changeCost);
4680                              //computeObjectiveValue();
4681                              createRim4(false);
4682                              // make sure duals are current
4683                              computeDuals(givenDuals);
4684                              checkDualSolution();
4685                              progress_.modifyObjective(-COIN_DBL_MAX);
4686#define DUAL_TRY_FASTER
4687#ifdef DUAL_TRY_FASTER
4688                              if (numberDualInfeasibilities_) {
4689#endif
4690                                   numberChanged_ = 1; // force something to happen
4691                                   lastCleaned = numberIterations_ - 1;
4692#ifdef DUAL_TRY_FASTER
4693                              } else {
4694                                   //double value = objectiveValue_;
4695                                   computeObjectiveValue(true);
4696                                   //printf("old %g new %g\n",value,objectiveValue_);
4697                                   //numberChanged_=1;
4698                              }
4699#endif
4700                         }
4701                         if (lastCleaned < numberIterations_ && numberTimesOptimal_ < 4 &&
4702                                   (numberChanged_ || (specialOptions_ & 4096) == 0)) {
4703#if CLP_CAN_HAVE_ZERO_OBJ
4704                           if ((specialOptions_&2097152)==0) {
4705#endif
4706                              doOriginalTolerance = 2;
4707                              numberTimesOptimal_++;
4708                              changeMade_++; // say something changed
4709                              if (numberTimesOptimal_ == 1) {
4710                                   dualTolerance_ = dblParam_[ClpDualTolerance];
4711                              } else {
4712                                   if (numberTimesOptimal_ == 2) {
4713                                        // better to have small tolerance even if slower
4714                                        factorization_->zeroTolerance(CoinMin(factorization_->zeroTolerance(), 1.0e-15));
4715                                   }
4716                                   dualTolerance_ = dblParam_[ClpDualTolerance];
4717                                   dualTolerance_ *= pow(2.0, numberTimesOptimal_ - 1);
4718                              }
4719                              cleanDuals = 2; // If nothing changed optimal else primal
4720#if CLP_CAN_HAVE_ZERO_OBJ
4721                           } else {
4722                             // no cost - skip checks
4723                             problemStatus_=0;
4724                           }
4725#endif
4726                         } else {
4727                              problemStatus_ = 0; // optimal
4728                              if (lastCleaned < numberIterations_ && numberChanged_) {
4729                                   handler_->message(CLP_SIMPLEX_GIVINGUP, messages_)
4730                                             << CoinMessageEol;
4731                              }
4732                         }
4733                    } else {
4734                         cleanDuals = 1;
4735                         if (doOriginalTolerance == 1) {
4736                              // check unbounded
4737                              // find a variable with bad dj
4738                              int iSequence;
4739                              int iChosen = -1;
4740                              if (!inCbc) {
4741                                   double largest = 100.0 * primalTolerance_;
4742                                   for (iSequence = 0; iSequence < numberRows_ + numberColumns_;
4743                                             iSequence++) {
4744                                        double djValue = dj_[iSequence];
4745                                        double originalLo = originalLower(iSequence);
4746                                        double originalUp = originalUpper(iSequence);
4747                                        if (fabs(djValue) > fabs(largest)) {
4748                                             if (getStatus(iSequence) != basic) {
4749                                                  if (djValue > 0 && originalLo < -1.0e20) {
4750                                                       if (djValue > fabs(largest)) {
4751                                                            largest = djValue;
4752                                                            iChosen = iSequence;
4753                                                       }
4754                                                  } else if (djValue < 0 && originalUp > 1.0e20) {
4755                                                       if (-djValue > fabs(largest)) {
4756                                                            largest = djValue;
4757                                                            iChosen = iSequence;
4758                                                       }
4759                                                  }
4760                                             }
4761                                        }
4762                                   }
4763                              }
4764                              if (iChosen >= 0) {
4765                                   int iSave = sequenceIn_;
4766                                   sequenceIn_ = iChosen;
4767                                   unpack(rowArray_[1]);
4768                                   sequenceIn_ = iSave;
4769                                   // if dual infeasibilities then must be free vector so add in dual
4770                                   if (numberDualInfeasibilities_) {
4771                                        if (fabs(changeCost) > 1.0e-5)
4772                                          COIN_DETAIL_PRINT(printf("Odd free/unbounded combo\n"));
4773                                        changeCost += cost_[iChosen];
4774                                   }
4775                                   problemStatus_ = checkUnbounded(rowArray_[1], rowArray_[0],
4776                                                                   changeCost);
4777                                   rowArray_[1]->clear();
4778                              } else {
4779                                   problemStatus_ = -3;
4780                              }
4781                              if (problemStatus_ == 2 && perturbation_ == 101) {
4782                                   perturbation_ = 102; // stop any perturbations
4783                                   cleanDuals = 1;
4784                                   createRim4(false);
4785                                   progress_.modifyObjective(-COIN_DBL_MAX);
4786                                   problemStatus_ = -1;
4787                              }
4788                              if (problemStatus_ == 2) {
4789                                   // it is unbounded - restore solution
4790                                   // but first add in changes to non-basic
4791                                   int iColumn;
4792                                   double * original = columnArray_[0]->denseVector();
4793                                   for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
4794                                        if(getColumnStatus(iColumn) != basic)
4795                                             ray_[iColumn] +=
4796                                                  saveColumnSolution[iColumn] - original[iColumn];
4797                                        columnActivityWork_[iColumn] = original[iColumn];
4798                                   }
4799                                   CoinMemcpyN(saveRowSolution, numberRows_,
4800                                               rowActivityWork_);
4801                              }
4802                         } else {
4803                              doOriginalTolerance = 2;
4804                              rowArray_[0]->clear();
4805                         }
4806                    }
4807                    delete [] saveColumnSolution;
4808                    delete [] saveRowSolution;
4809               }
4810               if (problemStatus_ == -4 || problemStatus_ == -5) {
4811                    // may be infeasible - or may be bounds are wrong
4812                    numberChangedBounds = changeBounds(0, NULL, changeCost);
4813                    needCleanFake = true;
4814                    /* Should this be here as makes no difference to being feasible.
4815                       But seems to make a difference to run times. */
4816                    if (perturbation_ == 101 && 0) {
4817                         perturbation_ = 102; // stop any perturbations
4818                         cleanDuals = 1;
4819                         numberChangedBounds = 1;
4820                         // make sure fake bounds are back
4821                         changeBounds(1, NULL, changeCost);
4822                         needCleanFake = true;
4823                         createRim4(false);
4824                         progress_.modifyObjective(-COIN_DBL_MAX);
4825                    }
4826                    if ((numberChangedBounds <= 0 || dualBound_ > 1.0e20 ||
4827                              (largestPrimalError_ > 1.0 && dualBound_ > 1.0e17)) &&
4828                              (numberPivots < 4 || sumPrimalInfeasibilities_ > 1.0e-6)) {
4829                         problemStatus_ = 1; // infeasible
4830                         if (perturbation_ == 101) {
4831                              perturbation_ = 102; // stop any perturbations
4832                              //cleanDuals=1;
4833                              //numberChangedBounds=1;
4834                              //createRim4(false);
4835                         }
4836                    } else {
4837                         problemStatus_ = -1; //iterate
4838                         cleanDuals = 1;
4839                         if (numberChangedBounds <= 0)
4840                              doOriginalTolerance = 2;
4841                         // and delete ray which has been created
4842                         delete [] ray_;
4843                         ray_ = NULL;
4844                    }
4845
4846               }
4847          } else {
4848               cleanDuals = 1;
4849          }
4850          if (problemStatus_ < 0) {
4851               if (doOriginalTolerance == 2) {
4852                    // put back original tolerance
4853                    lastCleaned = numberIterations_;
4854                    numberChanged_ = 0; // Number of variables with changed costs
4855                    handler_->message(CLP_DUAL_ORIGINAL, messages_)
4856                              << CoinMessageEol;
4857                    perturbation_ = 102; // stop any perturbations
4858#if 0
4859                    double * xcost = new double[numberRows_+numberColumns_];
4860                    double * xlower = new double[numberRows_+numberColumns_];
4861                    double * xupper = new double[numberRows_+numberColumns_];
4862                    double * xdj = new double[numberRows_+numberColumns_];
4863                    double * xsolution = new double[numberRows_+numberColumns_];
4864                    CoinMemcpyN(cost_, (numberRows_ + numberColumns_), xcost);
4865                    CoinMemcpyN(lower_, (numberRows_ + numberColumns_), xlower);
4866                    CoinMemcpyN(upper_, (numberRows_ + numberColumns_), xupper);
4867                    CoinMemcpyN(dj_, (numberRows_ + numberColumns_), xdj);
4868                    CoinMemcpyN(solution_, (numberRows_ + numberColumns_), xsolution);
4869#endif
4870                    createRim4(false);
4871                    progress_.modifyObjective(-COIN_DBL_MAX);
4872                    // make sure duals are current
4873                    computeDuals(givenDuals);
4874                    checkDualSolution();
4875#if 0
4876                    int i;
4877                    for (i = 0; i < numberRows_ + numberColumns_; i++) {
4878                         if (cost_[i] != xcost[i])
4879                              printf("** %d old cost %g new %g sol %g\n",
4880                                     i, xcost[i], cost_[i], solution_[i]);
4881                         if (lower_[i] != xlower[i])
4882                              printf("** %d old lower %g new %g sol %g\n",
4883                                     i, xlower[i], lower_[i], solution_[i]);
4884                         if (upper_[i] != xupper[i])
4885                              printf("** %d old upper %g new %g sol %g\n",
4886                                     i, xupper[i], upper_[i], solution_[i]);
4887                         if (dj_[i] != xdj[i])
4888                              printf("** %d old dj %g new %g sol %g\n",
4889                                     i, xdj[i], dj_[i], solution_[i]);
4890                         if (solution_[i] != xsolution[i])
4891                              printf("** %d old solution %g new %g sol %g\n",
4892                                     i, xsolution[i], solution_[i], solution_[i]);
4893                    }
4894                    //delete [] xcost;
4895                    //delete [] xupper;
4896                    //delete [] xlower;
4897                    //delete [] xdj;
4898                    //delete [] xsolution;
4899#endif
4900                    // put back bounds as they were if was optimal
4901                    if (doOriginalTolerance == 2 && cleanDuals != 2) {
4902                         changeMade_++; // say something changed
4903                         /* We may have already changed some bounds in this function
4904                            so save numberFake_ and add in.
4905
4906                            Worst that can happen is that we waste a bit of time  - but it must be finite.
4907                         */
4908                         //int saveNumberFake = numberFake_;
4909                         //resetFakeBounds(-1);
4910                         changeBounds(3, NULL, changeCost);
4911                         needCleanFake = true;
4912                         //numberFake_ += saveNumberFake;
4913                         //resetFakeBounds(-1);
4914                         cleanDuals = 2;
4915                         //cleanDuals=1;
4916                    }
4917#if 0
4918                    //int i;
4919                    for (i = 0; i < numberRows_ + numberColumns_; i++) {
4920                         if (cost_[i] != xcost[i])
4921                              printf("** %d old cost %g new %g sol %g\n",
4922                                     i, xcost[i], cost_[i], solution_[i]);
4923                         if (lower_[i] != xlower[i])
4924                              printf("** %d old lower %g new %g sol %g\n",
4925                                     i, xlower[i], lower_[i], solution_[i]);
4926                         if (upper_[i] != xupper[i])
4927                              printf("** %d old upper %g new %g sol %g\n",
4928                                     i, xupper[i], upper_[i], solution_[i]);
4929                         if (dj_[i] != xdj[i])
4930                              printf("** %d old dj %g new %g sol %g\n",
4931                                     i, xdj[i], dj_[i], solution_[i]);
4932                         if (solution_[i] != xsolution[i])
4933                              printf("** %d old solution %g new %g sol %g\n",
4934                                     i, xsolution[i], solution_[i], solution_[i]);
4935                    }
4936                    delete [] xcost;
4937                    delete [] xupper;
4938                    delete [] xlower;
4939                    delete [] xdj;
4940                    delete [] xsolution;
4941#endif
4942               }
4943               if (cleanDuals == 1 || (cleanDuals == 2 && !numberDualInfeasibilities_)) {
4944                    // make sure dual feasible
4945                    // look at all rows and columns
4946                    rowArray_[0]->clear();
4947                    columnArray_[0]->clear();
4948                    double objectiveChange = 0.0;
4949                    double savePrimalInfeasibilities = sumPrimalInfeasibilities_;
4950                    if (!numberIterations_) {
4951                      int nTotal = numberRows_ + numberColumns_;
4952                      if (arraysNotCreated) {
4953                        // create save arrays
4954                        delete [] saveStatus_;
4955                        delete [] savedSolution_;
4956                        saveStatus_ = new unsigned char [nTotal];
4957                        savedSolution_ = new double [nTotal];
4958                        arraysNotCreated = false;
4959                      }
4960                      // save arrays
4961                      CoinMemcpyN(status_, nTotal, saveStatus_);
4962                      CoinMemcpyN(rowActivityWork_,
4963                                  numberRows_, savedSolution_ + numberColumns_);
4964                      CoinMemcpyN(columnActivityWork_, numberColumns_, savedSolution_);
4965                    }
4966#if 0
4967                    double * xcost = new double[numberRows_+numberColumns_];
4968                    double * xlower = new double[numberRows_+numberColumns_];
4969                    double * xupper = new double[numberRows_+numberColumns_];
4970                    double * xdj = new double[numberRows_+numberColumns_];
4971                    double * xsolution = new double[numberRows_+numberColumns_];
4972                    CoinMemcpyN(cost_, (numberRows_ + numberColumns_), xcost);
4973                    CoinMemcpyN(lower_, (numberRows_ + numberColumns_), xlower);
4974                    CoinMemcpyN(upper_, (numberRows_ + numberColumns_), xupper);
4975                    CoinMemcpyN(dj_, (numberRows_ + numberColumns_), xdj);
4976                    CoinMemcpyN(solution_, (numberRows_ + numberColumns_), xsolution);
4977#endif
4978                    if (givenDuals)
4979                         dualTolerance_ = 1.0e50;
4980#if CLP_CAN_HAVE_ZERO_OBJ>1
4981                    if ((specialOptions_&2097152)==0) {
4982#endif
4983                    updateDualsInDual(rowArray_[0], columnArray_[0], rowArray_[1],
4984                                      0.0, objectiveChange, true);
4985#if CLP_CAN_HAVE_ZERO_OBJ>1
4986                    } else {
4987                      rowArray_[0]->clear();
4988                      rowArray_[1]->clear();
4989                      columnArray_[0]->clear();
4990                    }
4991#endif
4992                    dualTolerance_ = saveTolerance;
4993#if 0
4994                    int i;
4995                    for (i = 0; i < numberRows_ + numberColumns_; i++) {
4996                         if (cost_[i] != xcost[i])
4997                              printf("** %d old cost %g new %g sol %g\n",
4998                                     i, xcost[i], cost_[i], solution_[i]);
4999                         if (lower_[i] != xlower[i])
5000                              printf("** %d old lower %g new %g sol %g\n",
5001                                     i, xlower[i], lower_[i], solution_[i]);
5002                         if (upper_[i] != xupper[i])
5003                              printf("** %d old upper %g new %g sol %g\n",
5004                                     i, xupper[i], upper_[i], solution_[i]);
5005                         if (dj_[i] != xdj[i])
5006                              printf("** %d old dj %g new %g sol %g\n",
5007                                     i, xdj[i], dj_[i], solution_[i]);
5008                         if (solution_[i] != xsolution[i])
5009                              printf("** %d old solution %g new %g sol %g\n",
5010                                     i, xsolution[i], solution_[i], solution_[i]);
5011                    }
5012                    delete [] xcost;
5013                    delete [] xupper;
5014                    delete [] xlower;
5015                    delete [] xdj;
5016                    delete [] xsolution;
5017#endif
5018                    // for now - recompute all
5019                    gutsOfSolution(NULL, NULL);
5020                    if (givenDuals)
5021                         dualTolerance_ = 1.0e50;
5022#if CLP_CAN_HAVE_ZERO_OBJ>1
5023                    if ((specialOptions_&2097152)==0) {
5024#endif
5025                    updateDualsInDual(rowArray_[0], columnArray_[0], rowArray_[1],
5026                                      0.0, objectiveChange, true);
5027#if CLP_CAN_HAVE_ZERO_OBJ>1
5028                    } else {
5029                      rowArray_[0]->clear();
5030                      rowArray_[1]->clear();
5031                      columnArray_[0]->clear();
5032                    }
5033#endif
5034                    dualTolerance_ = saveTolerance;
5035                    if (!numberIterations_ && sumPrimalInfeasibilities_ >
5036                        1.0e5*(savePrimalInfeasibilities+1.0e3) &&
5037                        (moreSpecialOptions_ & (256|8192)) == 0) {
5038                      // Use primal
5039                      int nTotal = numberRows_ + numberColumns_;
5040                      CoinMemcpyN(saveStatus_, nTotal, status_);
5041                      CoinMemcpyN(savedSolution_ + numberColumns_ ,
5042                                  numberRows_, rowActivityWork_);
5043                      CoinMemcpyN(savedSolution_ ,
5044                                  numberColumns_, columnActivityWork_);
5045                      problemStatus_ = 10;
5046                      situationChanged = 0;
5047                    }
5048                    //assert(numberDualInfeasibilitiesWithoutFree_==0);
5049                    if (numberDualInfeasibilities_) {
5050                        if ((numberPrimalInfeasibilities_ || numberPivots)
5051                            && problemStatus_!=10) {
5052                              problemStatus_ = -1; // carry on as normal
5053                         } else {
5054                              problemStatus_ = 10; // try primal
5055#if COIN_DEVELOP>1
5056                              printf("returning at %d\n", __LINE__);
5057#endif
5058                         }
5059                    } else if (situationChanged == 2) {
5060                         problemStatus_ = -1; // carry on as normal
5061                         // need to reset bounds
5062                         changeBounds(3, NULL, changeCost);
5063                    }
5064                    situationChanged = 0;
5065               } else {
5066                    // iterate
5067                    if (cleanDuals != 2) {
5068                         problemStatus_ = -1;
5069                    } else {
5070                         problemStatus_ = 10; // try primal
5071#if COIN_DEVELOP>2
5072                         printf("returning at %d\n", __LINE__);
5073#endif
5074                    }
5075               }
5076          }
5077     }
5078     // unflag all variables (we may want to wait a bit?)
5079     if ((tentativeStatus != -2 && tentativeStatus != -1) && unflagVariables) {
5080          int iRow;
5081          int numberFlagged = 0;
5082          for (iRow = 0; iRow < numberRows_; iRow++) {
5083               int iPivot = pivotVariable_[iRow];
5084               if (flagged(iPivot)) {
5085                    numberFlagged++;
5086                    clearFlagged(iPivot);
5087               }
5088          }
5089#ifdef COIN_DEVELOP
5090          if (numberFlagged) {
5091               printf("unflagging %d variables - tentativeStatus %d probStat %d ninf %d nopt %d\n", numberFlagged, tentativeStatus,
5092                      problemStatus_, numberPrimalInfeasibilities_,
5093                      numberTimesOptimal_);
5094          }
5095#endif
5096          unflagVariables = numberFlagged > 0;
5097          if (numberFlagged && !numberPivots) {
5098               /* looks like trouble as we have not done any iterations.
5099               Try changing pivot tolerance then give it a few goes and give up */
5100               if (factorization_->pivotTolerance() < 0.9) {
5101                    factorization_->pivotTolerance(0.99);
5102                    problemStatus_ = -1;
5103               } else if (numberTimesOptimal_ < 3) {
5104                    numberTimesOptimal_++;
5105                    problemStatus_ = -1;
5106               } else {
5107                    unflagVariables = false;
5108                    //secondaryStatus_ = 1; // and say probably infeasible
5109                    if ((moreSpecialOptions_ & 256) == 0) {
5110                         // try primal
5111                         problemStatus_ = 10;
5112                    } else {
5113                         // almost certainly infeasible
5114                         problemStatus_ = 1;
5115                    }
5116#if COIN_DEVELOP>1
5117                    printf("returning at %d\n", __LINE__);
5118#endif
5119               }
5120          }
5121     }
5122     if (problemStatus_ < 0) {
5123          if (needCleanFake) {
5124               double dummyChangeCost = 0.0;
5125               changeBounds(3, NULL, dummyChangeCost);
5126          }
5127#if 0
5128          if (objectiveValue_ < lastObjectiveValue_ - 1.0e-8 *
5129                    CoinMax(fabs(objectivevalue_), fabs(lastObjectiveValue_))) {
5130          } else {
5131               lastObjectiveValue_ = objectiveValue_;
5132          }
5133#endif
5134          if (type == 0 || type == 1) {
5135               if (!type && arraysNotCreated) {
5136                    // create save arrays
5137                    delete [] saveStatus_;
5138                    delete [] savedSolution_;
5139                    saveStatus_ = new unsigned char [numberRows_+numberColumns_];
5140                    savedSolution_ = new double [numberRows_+numberColumns_];
5141               }
5142               // save arrays
5143               CoinMemcpyN(status_, numberColumns_ + numberRows_, saveStatus_);
5144               CoinMemcpyN(rowActivityWork_,
5145                           numberRows_, savedSolution_ + numberColumns_);
5146               CoinMemcpyN(columnActivityWork_, numberColumns_, savedSolution_);
5147               // save extra stuff
5148               int dummy;
5149               matrix_->generalExpanded(this, 5, dummy);
5150          }
5151          if (weightsSaved) {
5152               // restore weights (if saved) - also recompute infeasibility list
5153               if (!reallyBadProblems && (largestPrimalError_ < 100.0 || numberPivots > 10)) {
5154                    if (tentativeStatus > -3)
5155                         dualRowPivot_->saveWeights(this, (type < 2) ? 2 : 4);
5156                    else
5157                         dualRowPivot_->saveWeights(this, 3);
5158               } else {
5159                    // reset weights or scale back
5160                    dualRowPivot_->saveWeights(this, 6);
5161               }
5162          } else if (weightsSaved2 && numberPrimalInfeasibilities_) {
5163               dualRowPivot_->saveWeights(this, 3);
5164          }
5165     }
5166     // see if cutoff reached
5167     double limit = 0.0;
5168     getDblParam(ClpDualObjectiveLimit, limit);
5169#if 0
5170     if(fabs(limit) < 1.0e30 && objectiveValue()*optimizationDirection_ >
5171               limit + 100.0) {
5172          printf("lim %g obj %g %g - wo perturb %g sum dual %g\n",
5173                 limit, objectiveValue_, objectiveValue(), computeInternalObjectiveValue(), sumDualInfeasibilities_);
5174     }
5175#endif
5176     if(fabs(limit) < 1.0e30 && objectiveValue()*optimizationDirection_ >
5177               limit && !numberAtFakeBound()) {
5178          bool looksInfeasible = !numberDualInfeasibilities_;
5179          if (objectiveValue()*optimizationDirection_ > limit + fabs(0.1 * limit) + 1.0e2 * sumDualInfeasibilities_ + 1.0e4 &&
5180                    sumDualInfeasibilities_ < largestDualError_ && numberIterations_ > 0.5 * numberRows_ + 1000)
5181               looksInfeasible = true;
5182          if (looksInfeasible) {
5183               // Even if not perturbed internal costs may have changed
5184               // be careful
5185               if (true || numberIterations_) {
5186                    if(computeInternalObjectiveValue() > limit) {
5187                         problemStatus_ = 1;
5188                         secondaryStatus_ = 1; // and say was on cutoff
5189                    }
5190               } else {
5191                    problemStatus_ = 1;
5192                    secondaryStatus_ = 1; // and say was on cutoff
5193               }
5194          }
5195     }
5196     // If we are in trouble and in branch and bound give up
5197     if ((specialOptions_ & 1024) != 0) {
5198          int looksBad = 0;
5199          if (largestPrimalError_ * largestDualError_ > 1.0e2) {
5200               looksBad = 1;
5201          } else if (largestPrimalError_ > 1.0e-2
5202                     && objectiveValue_ > CoinMin(1.0e15, 1.0e3 * limit)) {
5203               looksBad = 2;
5204          }
5205          if (looksBad) {
5206               if (factorization_->pivotTolerance() < 0.9) {
5207                    // up tolerance
5208                    factorization_->pivotTolerance(CoinMin(factorization_->pivotTolerance() * 1.05 + 0.02, 0.91));
5209               } else if (numberIterations_ > 10000) {
5210                    if (handler_->logLevel() > 2)
5211                         printf("bad dual - saying infeasible %d\n", looksBad);
5212                    problemStatus_ = 1;
5213                    secondaryStatus_ = 1; // and say was on cutoff
5214               } else if (largestPrimalError_ > 1.0e5) {
5215                    {
5216                      //int iBigB = -1;
5217                         double bigB = 0.0;
5218                         //int iBigN = -1;
5219                         double bigN = 0.0;
5220                         for (int i = 0; i < numberRows_ + numberColumns_; i++) {
5221                              double value = fabs(solution_[i]);
5222                              if (getStatus(i) == basic) {
5223                                   if (value > bigB) {
5224                                        bigB = value;
5225                                        //iBigB = i;
5226                                   }
5227                              } else {
5228                                   if (value > bigN) {
5229                                        bigN = value;
5230                                        //iBigN = i;
5231                                   }
5232                              }
5233                         }
5234#ifdef CLP_INVESTIGATE
5235                         if (bigB > 1.0e8 || bigN > 1.0e8) {
5236                              if (handler_->logLevel() > 0)
5237                                   printf("it %d - basic %d %g, nonbasic %d %g\n",
5238                                          numberIterations_, iBigB, bigB, iBigN, bigN);
5239                         }
5240#endif
5241                    }
5242#if COIN_DEVELOP!=2
5243                    if (handler_->logLevel() > 2)
5244#endif
5245                         printf("bad dual - going to primal %d %g\n", looksBad, largestPrimalError_);
5246                    allSlackBasis(true);
5247                    problemStatus_ = 10;
5248               }
5249          }
5250     }
5251     if (problemStatus_ < 0 && !changeMade_) {
5252          problemStatus_ = 4; // unknown
5253     }
5254     lastGoodIteration_ = numberIterations_;
5255     if (numberIterations_ > lastBadIteration_ + 100)
5256          moreSpecialOptions_ &= ~16; // clear check accuracy flag
5257     if (problemStatus_ < 0) {
5258          sumDualInfe