source: trunk/Clp/src/ClpSimplexDual.cpp

Last change on this file was 2544, checked in by unxusr, 4 weeks ago

interval between messages in Simplex implementations

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