source: stable/1.15/Clp/src/ClpSimplexDual.cpp @ 1989

Last change on this file since 1989 was 1931, checked in by stefan, 7 years ago

sync with trunk rev 1930

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