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

Last change on this file since 1706 was 1706, checked in by forrest, 9 years ago

fix for number iterations in ClpSolve? and check if gone feasible for extreme cases

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