source: stable/1.16/Clp/src/ClpSimplexDual.cpp @ 2205

Last change on this file since 2205 was 2205, checked in by forrest, 4 years ago

move delete in Idiot and try and improve free variables in dual

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