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

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

minor changes to implement Aboca

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