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

Last change on this file since 2385 was 2385, checked in by unxusr, 4 months ago

formatting

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