# source:stable/1.16/Clp/src/ClpSimplexDual.cpp@2102

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

try and trap problem going awry

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