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

Last change on this file since 2440 was 2440, checked in by forrest, 7 months ago

need changes in svn as well as git

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 295.8 KB
Line 
1/* $Id: ClpSimplexDual.cpp 2440 2019-04-04 14:44:30Z forrest $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6/* Notes on implementation of dual simplex algorithm.
7
8   When dual feasible:
9
10   If primal feasible, we are optimal.  Otherwise choose an infeasible
11   basic variable to leave basis (normally going to nearest bound) (B).  We
12   now need to find an incoming variable which will leave problem
13   dual feasible so we get the row of the tableau corresponding to
14   the basic variable (with the correct sign depending if basic variable
15   above or below feasibility region - as that affects whether reduced
16   cost on outgoing variable has to be positive or negative).
17
18   We now perform a ratio test to determine which incoming variable will
19   preserve dual feasibility (C).  If no variable found then problem
20   is infeasible (in primal sense).  If there is a variable, we then
21   perform pivot and repeat.  Trivial?
22
23   -------------------------------------------
24
25   A) How do we get dual feasible?  If all variables have bounds then
26   it is trivial to get feasible by putting non-basic variables to
27   correct bounds.  OSL did not have a phase 1/phase 2 approach but
28   instead effectively put fake bounds on variables and this is the
29   approach here, although I had hoped to make it cleaner.
30
31   If there is a weight of X on getting dual feasible:
32     Non-basic variables with negative reduced costs are put to
33     lesser of their upper bound and their lower bound + X.
34     Similarly, mutatis mutandis, for positive reduced costs.
35
36   Free variables should normally be in basis, otherwise I have
37   coding which may be able to come out (and may not be correct).
38
39   In OSL, this weight was changed heuristically, here at present
40   it is only increased if problem looks finished.  If problem is
41   feasible I check for unboundedness.  If not unbounded we
42   could play with going into primal.  As long as weights increase
43   any algorithm would be finite.
44
45   B) Which outgoing variable to choose is a virtual base class.
46   For difficult problems steepest edge is preferred while for
47   very easy (large) problems we will need partial scan.
48
49   C) Sounds easy, but this is hardest part of algorithm.
50      1) Instead of stopping at first choice, we may be able
51      to flip that variable to other bound and if objective
52      still improving choose again.  These mini iterations can
53      increase speed by orders of magnitude but we may need to
54      go to more of a bucket choice of variable rather than looking
55      at them one by one (for speed).
56      2) Accuracy.  Reduced costs may be of wrong sign but less than
57      tolerance.  Pivoting on these makes objective go backwards.
58      OSL modified cost so a zero move was made, Gill et al
59      (in primal analogue) modified so a strictly positive move was
60      made.  It is not quite as neat in dual but that is what we
61      try and do.  The two problems are that re-factorizations can
62      change reduced costs above and below tolerances and that when
63      finished we need to reset costs and try again.
64      3) Degeneracy.  Gill et al helps but may not be enough.  We
65      may need more.  Also it can improve speed a lot if we perturb
66      the costs significantly.
67
68  References:
69     Forrest and Goldfarb, Steepest-edge simplex algorithms for
70       linear programming - Mathematical Programming 1992
71     Forrest and Tomlin, Implementing the simplex method for
72       the Optimization Subroutine Library - IBM Systems Journal 1992
73     Gill, Murray, Saunders, Wright A Practical Anti-Cycling
74       Procedure for Linear and Nonlinear Programming SOL report 1988
75
76
77  TODO:
78
79  a) Better recovery procedures.  At present I never check on forward
80     progress.  There is checkpoint/restart with reducing
81     re-factorization frequency, but this is only on singular
82     factorizations.
83  b) Fast methods for large easy problems (and also the option for
84     the code to automatically choose which method).
85  c) We need to be able to stop in various ways for OSI - this
86     is fairly easy.
87
88 */
89#ifdef COIN_DEVELOP
90#undef COIN_DEVELOP
91#define COIN_DEVELOP 2
92#endif
93
94#include "CoinPragma.hpp"
95
96#include <math.h>
97
98#include "CoinHelperFunctions.hpp"
99#include "ClpHelperFunctions.hpp"
100#if ABOCA_LITE
101// 2 is owner of abcState_
102#define ABCSTATE_LITE 2
103#endif
104//#define FAKE_CILK
105#include "ClpSimplexDual.hpp"
106#include "ClpEventHandler.hpp"
107#include "ClpFactorization.hpp"
108#include "CoinPackedMatrix.hpp"
109#include "CoinIndexedVector.hpp"
110#include "CoinFloatEqual.hpp"
111#include "ClpDualRowDantzig.hpp"
112#include "ClpMessage.hpp"
113#include "ClpLinearObjective.hpp"
114#include <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                      (moreSpecialOptions_ & 256) == 0)
1947                    objective_ = new ClpLinearObjective(NULL, numberColumns_);
1948                }
1949              }
1950#else
1951              if (sumPrimalInfeasibilities_ < 1.0e-3 || sumDualInfeasibilities_ > 1.0e-6) {
1952#ifdef COIN_DEVELOP
1953                printf("at %d - primal %d %g - dual %d %g fake %d weight %g - pivs %d\n",
1954                  __LINE__, numberPrimalInfeasibilities_,
1955                  sumPrimalInfeasibilities_,
1956                  numberDualInfeasibilities_, sumDualInfeasibilities_,
1957                  numberFake_, dualBound_, factorization_->pivots());
1958#endif
1959                if ((specialOptions_ & 1024) != 0 && factorization_->pivots()) {
1960                  problemStatus_ = 10;
1961#if COIN_DEVELOP > 1
1962                  printf("returning at %d\n", __LINE__);
1963#endif
1964                  // Get rid of objective
1965                  if ((specialOptions_ & 16384) == 0 &&
1966                      (moreSpecialOptions_ & 256) == 0)
1967                    objective_ = new ClpLinearObjective(NULL, numberColumns_);
1968                }
1969              }
1970#endif
1971              rowArray_[0]->clear();
1972              columnArray_[0]->clear();
1973              returnCode = 1;
1974              break;
1975            }
1976          }
1977          // If special option set - put off as long as possible
1978          if ((specialOptions_ & 64) == 0 || (moreSpecialOptions_ & 64) != 0) {
1979            if (factorization_->pivots() == 0)
1980              problemStatus_ = -4; //say looks infeasible
1981          } else {
1982            // flag
1983            char x = isColumn(sequenceOut_) ? 'C' : 'R';
1984            handler_->message(CLP_SIMPLEX_FLAG, messages_)
1985              << x << sequenceWithin(sequenceOut_)
1986              << CoinMessageEol;
1987#ifdef COIN_DEVELOP
1988            printf("flag c\n");
1989#endif
1990            setFlagged(sequenceOut_);
1991            if (!factorization_->pivots()) {
1992              rowArray_[0]->clear();
1993              columnArray_[0]->clear();
1994              continue;
1995            }
1996          }
1997        }
1998        acceptablePivot_ = fabs(acceptablePivot_);
1999        if (factorization_->pivots() < 5 && acceptablePivot_ > 1.0e-8)
2000          acceptablePivot_ = 1.0e-8;
2001        rowArray_[0]->clear();
2002        columnArray_[0]->clear();
2003        returnCode = 1;
2004        break;
2005      }
2006    } else {
2007#ifdef CLP_INVESTIGATE_SERIAL
2008      z_thinks = 0;
2009#endif
2010      // no pivot row
2011#ifdef CLP_DEBUG
2012      if (handler_->logLevel() & 32)
2013        printf("** no row pivot\n");
2014#endif
2015      // If in branch and bound try and get rid of fixed variables
2016      if ((specialOptions_ & 1024) != 0 && CLEAN_FIXED) {
2017        assert(!candidateList);
2018        candidateList = new int[numberRows_];
2019        int iRow;
2020        for (iRow = 0; iRow < numberRows_; iRow++) {
2021          int iPivot = pivotVariable_[iRow];
2022          if (flagged(iPivot) || !pivoted(iPivot))
2023            continue;
2024          assert(iPivot < numberColumns_ && lower_[iPivot] == upper_[iPivot]);
2025          candidateList[numberCandidates++] = iRow;
2026        }
2027        // and set first candidate
2028        if (!numberCandidates) {
2029          delete[] candidateList;
2030          candidateList = NULL;
2031          int iRow;
2032          for (iRow = 0; iRow < numberRows_; iRow++) {
2033            int iPivot = pivotVariable_[iRow];
2034            clearPivoted(iPivot);
2035          }
2036        } else {
2037          ifValuesPass = 2;
2038          continue;
2039        }
2040      }
2041      int numberPivots = factorization_->pivots();
2042      bool specialCase;
2043      int useNumberFake;
2044      returnCode = 0;
2045      if (numberPivots <= CoinMax(dontFactorizePivots_, 20) && (specialOptions_ & 2048) != 0 && (true || !numberChanged_ || perturbation_ == 101)
2046        && dualBound_ >= 1.0e8) {
2047        specialCase = true;
2048        // as dual bound high - should be okay
2049        useNumberFake = 0;
2050      } else {
2051        specialCase = false;
2052        useNumberFake = numberFake_;
2053      }
2054      if (!numberPivots || specialCase) {
2055        if (numberPrimalInfeasibilities_ && problemStatus_ == -1)
2056          problemStatus_ = -4;
2057        // may have crept through - so may be optimal
2058        // check any flagged variables
2059        int iRow;
2060        for (iRow = 0; iRow < numberRows_; iRow++) {
2061          int iPivot = pivotVariable_[iRow];
2062          if (flagged(iPivot))
2063            break;
2064        }
2065        if (iRow < numberRows_ && numberPivots) {
2066          // try factorization
2067          returnCode = -2;
2068        }
2069
2070        if (useNumberFake || numberDualInfeasibilities_) {
2071          // may be dual infeasible
2072          if ((specialOptions_ & 1024) == 0)
2073            problemStatus_ = -5;
2074          else if (!useNumberFake && numberPrimalInfeasibilities_
2075            && !numberPivots)
2076            problemStatus_ = 1;
2077        } else {
2078          if (iRow < numberRows_) {
2079#ifdef COIN_DEVELOP
2080            std::cout << "Flagged variables at end - infeasible?" << std::endl;
2081            printf("Probably infeasible - pivot was %g\n", alpha_);
2082#endif
2083            //if (fabs(alpha_)<1.0e-4) {
2084            //problemStatus_=1;
2085            //} else {
2086#ifdef CLP_DEBUG
2087            abort();
2088#endif
2089            //}
2090            problemStatus_ = -5;
2091          } else {
2092            problemStatus_ = 0;
2093#ifndef CLP_CHECK_NUMBER_PIVOTS
2094#define CLP_CHECK_NUMBER_PIVOTS 10
2095#endif
2096#if CLP_CHECK_NUMBER_PIVOTS < 20
2097            if (numberPivots > CLP_CHECK_NUMBER_PIVOTS) {
2098#ifndef NDEBUG_CLP
2099              int nTotal = numberRows_ + numberColumns_;
2100              double *comp = CoinCopyOfArray(solution_, nTotal);
2101#endif
2102              computePrimals(rowActivityWork_, columnActivityWork_);
2103#ifndef NDEBUG_CLP
2104              double largest = 1.0e-5;
2105              int bad = -1;
2106              for (int i = 0; i < nTotal; i++) {
2107                double value = solution_[i];
2108                double larger = CoinMax(fabs(value), fabs(comp[i]));
2109                double tol = 1.0e-5 + 1.0e-5 * larger;
2110                double diff = fabs(value - comp[i]);
2111                if (diff - tol > largest) {
2112                  bad = i;
2113                  largest = diff - tol;
2114                }
2115              }
2116              if (bad >= 0)
2117                COIN_DETAIL_PRINT(printf("bad %d old %g new %g\n", bad, comp[bad], solution_[bad]));
2118#endif
2119              checkPrimalSolution(rowActivityWork_, columnActivityWork_);
2120              if (numberPrimalInfeasibilities_) {
2121#ifdef CLP_INVESTIGATE
2122                printf("XXX Infeas ? %d inf summing to %g\n", numberPrimalInfeasibilities_,
2123                  sumPrimalInfeasibilities_);
2124#endif
2125                problemStatus_ = -1;
2126                returnCode = -2;
2127              }
2128#ifndef NDEBUG_CLP
2129              memcpy(solution_, comp, nTotal * sizeof(double));
2130              delete[] comp;
2131#endif
2132            }
2133#endif
2134            if (!problemStatus_) {
2135              // make it look OK
2136              numberPrimalInfeasibilities_ = 0;
2137              sumPrimalInfeasibilities_ = 0.0;
2138              numberDualInfeasibilities_ = 0;
2139              sumDualInfeasibilities_ = 0.0;
2140              // May be perturbed
2141              if (perturbation_ == 101 || numberChanged_) {
2142                numberChanged_ = 0; // Number of variables with changed costs
2143                perturbation_ = 102; // stop any perturbations
2144                //double changeCost;
2145                //changeBounds(1,NULL,changeCost);
2146                createRim4(false);
2147                // make sure duals are current
2148                computeDuals(givenDuals);
2149                checkDualSolution();
2150                progress_.modifyObjective(-COIN_DBL_MAX);
2151                if (numberDualInfeasibilities_) {
2152                  problemStatus_ = 10; // was -3;
2153                } else {
2154                  computeObjectiveValue(true);
2155                }
2156              } else if (numberPivots) {
2157                computeObjectiveValue(true);
2158              }
2159              if (numberPivots < -1000) {
2160                // objective may be wrong
2161                objectiveValue_ = innerProduct(cost_, numberColumns_ + numberRows_, solution_);
2162                objectiveValue_ += objective_->nonlinearOffset();
2163                objectiveValue_ /= (objectiveScale_ * rhsScale_);
2164                if ((specialOptions_ & 16384) == 0) {
2165                  // and dual_ may be wrong (i.e. for fixed or basic)
2166                  CoinIndexedVector *arrayVector = rowArray_[1];
2167                  arrayVector->clear();
2168                  int iRow;
2169                  double *array = arrayVector->denseVector();
2170                  /* Use dual_ instead of array
2171                                                Even though dual_ is only numberRows_ long this is
2172                                                okay as gets permuted to longer rowArray_[2]
2173                                             */
2174                  arrayVector->setDenseVector(dual_);
2175                  int *index = arrayVector->getIndices();
2176                  int number = 0;
2177                  for (iRow = 0; iRow < numberRows_; iRow++) {
2178                    int iPivot = pivotVariable_[iRow];
2179                    double value = cost_[iPivot];
2180                    dual_[iRow] = value;
2181                    if (value) {
2182                      index[number++] = iRow;
2183                    }
2184                  }
2185                  arrayVector->setNumElements(number);
2186                  // Extended duals before "updateTranspose"
2187                  matrix_->dualExpanded(this, arrayVector, NULL, 0);
2188                  // Btran basic costs
2189                  rowArray_[2]->clear();
2190                  factorization_->updateColumnTranspose(rowArray_[2], arrayVector);
2191                  // and return vector
2192                  arrayVector->setDenseVector(array);
2193                }
2194              }
2195              sumPrimalInfeasibilities_ = 0.0;
2196            }
2197            if ((specialOptions_ & (1024 + 16384)) != 0 && !problemStatus_) {
2198              CoinIndexedVector *arrayVector = rowArray_[1];
2199              arrayVector->clear();
2200              double *rhs = arrayVector->denseVector();
2201              times(1.0, solution_, rhs);
2202#ifdef CHECK_ACCURACY
2203              bool bad = false;
2204#endif
2205              bool bad2 = false;
2206              int i;
2207              for (i = 0; i < numberRows_; i++) {
2208                if (rhs[i] < rowLowerWork_[i] - primalTolerance_ || rhs[i] > rowUpperWork_[i] + primalTolerance_) {
2209                  bad2 = true;
2210#ifdef CHECK_ACCURACY
2211                  printf("row %d out of bounds %g, %g correct %g bad %g\n", i,
2212                    rowLowerWork_[i], rowUpperWork_[i],
2213                    rhs[i], rowActivityWork_[i]);
2214#endif
2215                } else if (fabs(rhs[i] - rowActivityWork_[i]) > 1.0e-3) {
2216#ifdef CHECK_ACCURACY
2217                  bad = true;
2218                  printf("row %d correct %g bad %g\n", i, rhs[i], rowActivityWork_[i]);
2219#endif
2220                }
2221                rhs[i] = 0.0;
2222              }
2223              for (i = 0; i < numberColumns_; i++) {
2224                if (solution_[i] < columnLowerWork_[i] - primalTolerance_ || solution_[i] > columnUpperWork_[i] + primalTolerance_) {
2225                  bad2 = true;
2226#ifdef CHECK_ACCURACY
2227                  printf("column %d out of bounds %g, %g correct %g bad %g\n", i,
2228                    columnLowerWork_[i], columnUpperWork_[i],
2229                    solution_[i], columnActivityWork_[i]);
2230#endif
2231                }
2232              }
2233              if (bad2) {
2234                problemStatus_ = -3;
2235                returnCode = -2;
2236                // Force to re-factorize early next time
2237                int numberPivots = factorization_->pivots();
2238                forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
2239              }
2240            }
2241          }
2242        }
2243      } else {
2244        problemStatus_ = -3;
2245        returnCode = -2;
2246        // Force to re-factorize early next time
2247        int numberPivots = factorization_->pivots();
2248        forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
2249      }
2250      break;
2251    }
2252  }
2253  if (givenDuals) {
2254    CoinMemcpyN(dj_, numberRows_ + numberColumns_, givenDuals);
2255    // get rid of any values pass array
2256    delete[] candidateList;
2257  }
2258  delete[] dubiousWeights;
2259#ifdef CLP_REPORT_PROGRESS
2260  if (ixxxxxx > ixxyyyy - 5) {
2261    int nTotal = numberColumns_ + numberRows_;
2262    double oldObj = 0.0;
2263    double newObj = 0.0;
2264    for (int i = 0; i < nTotal; i++) {
2265      if (savePSol[i])
2266        oldObj += savePSol[i] * saveCost[i];
2267      if (solution_[i])
2268        newObj += solution_[i] * cost_[i];
2269      bool printIt = false;
2270      if (cost_[i] != saveCost[i])
2271        printIt = true;
2272      if (status_[i] != saveStat[i])
2273        printIt = true;
2274      if (printIt)
2275        printf("%d old %d cost %g sol %g, new %d cost %g sol %g\n",
2276          i, saveStat[i], saveCost[i], savePSol[i],
2277          status_[i], cost_[i], solution_[i]);
2278      // difference
2279      savePSol[i] = solution_[i] - savePSol[i];
2280    }
2281    printf("exit pivots %d, old obj %g new %g\n",
2282      factorization_->pivots(),
2283      oldObj, newObj);
2284    memset(saveDj, 0, numberRows_ * sizeof(double));
2285    times(1.0, savePSol, saveDj);
2286    double largest = 1.0e-6;
2287    int k = -1;
2288    for (int i = 0; i < numberRows_; i++) {
2289      saveDj[i] -= savePSol[i + numberColumns_];
2290      if (fabs(saveDj[i]) > largest) {
2291        largest = fabs(saveDj[i]);
2292        k = i;
2293      }
2294    }
2295    if (k >= 0)
2296      printf("Not null %d %g\n", k, largest);
2297  }
2298  delete[] savePSol;
2299  delete[] saveDj;
2300  delete[] saveCost;
2301  delete[] saveStat;
2302#endif
2303  return returnCode;
2304}
2305#if ABOCA_LITE
2306static void
2307updateDualBit(clpTempInfo &info)
2308{
2309  int numberInfeasibilities = 0;
2310  double tolerance = info.tolerance;
2311  double theta = info.theta;
2312  double *COIN_RESTRICT reducedCost = info.reducedCost;
2313  const double *COIN_RESTRICT lower = info.lower;
2314  const double *COIN_RESTRICT upper = info.upper;
2315  double *COIN_RESTRICT work = info.work;
2316  int number = info.numberToDo;
2317  int *COIN_RESTRICT which = info.which;
2318  const unsigned char *COIN_RESTRICT statusArray = info.status;
2319  double multiplier[] = { -1.0, 1.0 };
2320  for (int i = 0; i < number; i++) {
2321    int iSequence = which[i];
2322    double alphaI = work[i];
2323    work[i] = 0.0;
2324
2325    int iStatus = (statusArray[iSequence] & 3) - 1;
2326    if (iStatus) {
2327      double value = reducedCost[iSequence] - theta * alphaI;
2328      reducedCost[iSequence] = value;
2329      //printf("xx %d %.18g\n",iSequence,reducedCost[iSequence]);
2330      double mult = multiplier[iStatus - 1];
2331      value *= mult;
2332      // skip if free
2333      if (value < -tolerance && iStatus > 0) {
2334        // flipping bounds
2335        double movement = mult * (upper[iSequence] - lower[iSequence]);
2336        work[numberInfeasibilities] = movement;
2337        which[numberInfeasibilities++] = iSequence;
2338      }
2339    }
2340  }
2341  info.numberInfeasibilities = numberInfeasibilities;
2342}
2343#endif
2344/* The duals are updated by the given arrays.
2345   Returns number of infeasibilities.
2346   rowArray and columnarray will have flipped
2347   The output vector has movement (row length array) */
2348int ClpSimplexDual::updateDualsInDual(CoinIndexedVector *rowArray,
2349  CoinIndexedVector *columnArray,
2350  CoinIndexedVector *outputArray,
2351  double theta,
2352  double &objectiveChange,
2353  bool fullRecompute)
2354{
2355
2356  outputArray->clear();
2357
2358  int numberInfeasibilities = 0;
2359  int numberRowInfeasibilities = 0;
2360
2361  // get a tolerance
2362  double tolerance = dualTolerance_;
2363  // we can't really trust infeasibilities if there is dual error
2364  double error = CoinMin(1.0e-2, largestDualError_);
2365  // allow tolerance at least slightly bigger than standard
2366  tolerance = tolerance + error;
2367
2368  double changeObj = 0.0;
2369
2370  // Coding is very similar but we can save a bit by splitting
2371  // Do rows
2372  if (!fullRecompute) {
2373    int i;
2374    double *COIN_RESTRICT reducedCost = djRegion(0);
2375    const double *COIN_RESTRICT lower = lowerRegion(0);
2376    const double *COIN_RESTRICT upper = upperRegion(0);
2377    const double *COIN_RESTRICT cost = costRegion(0);
2378    double *COIN_RESTRICT work;
2379    int number;
2380    int *COIN_RESTRICT which;
2381    const unsigned char *COIN_RESTRICT statusArray = status_ + numberColumns_;
2382    assert(rowArray->packedMode());
2383    work = rowArray->denseVector();
2384    number = rowArray->getNumElements();
2385    which = rowArray->getIndices();
2386    double multiplier[] = { -1.0, 1.0 };
2387    for (i = 0; i < number; i++) {
2388      int iSequence = which[i];
2389      double alphaI = work[i];
2390      work[i] = 0.0;
2391      int iStatus = (statusArray[iSequence] & 3) - 1;
2392      if (iStatus) {
2393        double value = reducedCost[iSequence] - theta * alphaI;
2394        // NO - can have free assert (iStatus>0);
2395        reducedCost[iSequence] = value;
2396        double mult = multiplier[iStatus - 1];
2397        value *= mult;
2398        // skip if free
2399        if (value < -tolerance && iStatus > 0) {
2400          // flipping bounds
2401          double movement = mult * (lower[iSequence] - upper[iSequence]);
2402          which[numberInfeasibilities++] = iSequence;
2403#ifndef NDEBUG
2404          if (fabs(movement) >= 1.0e30)
2405            resetFakeBounds(-1000 - iSequence);
2406#endif
2407#ifdef CLP_DEBUG
2408          if ((handler_->logLevel() & 32))
2409            printf("%d %d, new dj %g, alpha %g, movement %g\n",
2410              0, iSequence, value, alphaI, movement);
2411#endif
2412          changeObj -= movement * cost[iSequence];
2413          outputArray->quickAdd(iSequence, movement);
2414        }
2415      }
2416    }
2417    // Do columns
2418    reducedCost = djRegion(1);
2419    lower = lowerRegion(1);
2420    upper = upperRegion(1);
2421    cost = costRegion(1);
2422    // set number of infeasibilities in row array
2423    numberRowInfeasibilities = numberInfeasibilities;
2424    rowArray->setNumElements(numberInfeasibilities);
2425    numberInfeasibilities = 0;
2426    work = columnArray->denseVector();
2427    number = columnArray->getNumElements();
2428    which = columnArray->getIndices();
2429    if ((moreSpecialOptions_ & 8) != 0) {
2430      const unsigned char *COIN_RESTRICT statusArray = status_;
2431#if ABOCA_LITE
2432      int numberThreads = abcState();
2433      if (numberThreads) {
2434        clpTempInfo info[ABOCA_LITE];
2435        int chunk = (number + numberThreads - 1) / numberThreads;
2436        int n = 0;
2437        int *whichX = which;
2438        for (i = 0; i < numberThreads; i++) {
2439          info[i].theta = theta;
2440          info[i].tolerance = tolerance;
2441          info[i].reducedCost = reducedCost;
2442          info[i].lower = lower;
2443          info[i].upper = upper;
2444          info[i].status = statusArray;
2445          info[i].which = which + n;
2446          info[i].work = work + n;
2447          info[i].numberToDo = CoinMin(chunk, number - n);
2448          n += chunk;
2449        }
2450        for (i = 0; i < numberThreads; i++) {
2451          cilk_spawn updateDualBit(info[i]);
2452        }
2453        cilk_sync;
2454        for (i = 0; i < numberThreads; i++) {
2455          int n = info[i].numberInfeasibilities;
2456          double *workV = info[i].work;
2457          int *whichV = info[i].which;
2458          for (int j = 0; j < n; j++) {
2459            int iSequence = whichV[j];
2460            double movement = workV[j];
2461            workV[j] = 0.0;
2462            whichX[numberInfeasibilities++] = iSequence;
2463#ifndef NDEBUG
2464            if (fabs(movement) >= 1.0e30)
2465              resetFakeBounds(-1000 - iSequence);
2466#endif
2467            changeObj += movement * cost[iSequence];
2468            matrix_->add(this, outputArray, iSequence, movement);
2469          }
2470        }
2471      } else {
2472#endif
2473        for (i = 0; i < number; i++) {
2474          int iSequence = which[i];
2475          double alphaI = work[i];
2476          work[i] = 0.0;
2477
2478          int iStatus = (statusArray[iSequence] & 3) - 1;
2479          if (iStatus) {
2480            double value = reducedCost[iSequence] - theta * alphaI;
2481            assert(iStatus > 0);
2482            reducedCost[iSequence] = value;
2483            //printf("xx %d %.18g\n",iSequence,reducedCost[iSequence]);
2484            double mult = multiplier[iStatus - 1];
2485            value *= mult;
2486            // skip if free
2487            if (value < -tolerance && iStatus > 0) {
2488              // flipping bounds
2489              double movement = mult * (upper[iSequence] - lower[iSequence]);
2490              which[numberInfeasibilities++] = iSequence;
2491#ifndef NDEBUG
2492              if (fabs(movement) >= 1.0e30)
2493                resetFakeBounds(-1000 - iSequence);
2494#endif
2495#ifdef CLP_DEBUG
2496              if ((handler_->logLevel() & 32))
2497                printf("%d %d, new dj %g, alpha %g, movement %g\n",
2498                  1, iSequence, value, alphaI, movement);
2499#endif
2500              changeObj += movement * cost[iSequence];
2501              matrix_->add(this, outputArray, iSequence, movement);
2502            }
2503          }
2504        }
2505#if ABOCA_LITE
2506      }
2507#endif
2508    } else {
2509      for (i = 0; i < number; i++) {
2510        int iSequence = which[i];
2511        double alphaI = work[i];
2512        work[i] = 0.0;
2513
2514        Status status = getStatus(iSequence);
2515        if (status == atLowerBound) {
2516          double value = reducedCost[iSequence] - theta * alphaI;
2517          reducedCost[iSequence] = value;
2518          double movement = 0.0;
2519
2520          if (value < -tolerance) {
2521            // to upper bound
2522            which[numberInfeasibilities++] = iSequence;
2523            movement = upper[iSequence] - lower[iSequence];
2524#ifndef NDEBUG
2525            if (fabs(movement) >= 1.0e30)
2526              resetFakeBounds(-1000 - iSequence);
2527#endif
2528#ifdef CLP_DEBUG
2529            if ((handler_->logLevel() & 32))
2530              printf("%d %d, new dj %g, alpha %g, movement %g\n",
2531                1, iSequence, value, alphaI, movement);
2532#endif
2533            changeObj += movement * cost[iSequence];
2534            matrix_->add(this, outputArray, iSequence, movement);
2535          }
2536        } else if (status == atUpperBound) {
2537          double value = reducedCost[iSequence] - theta * alphaI;
2538          reducedCost[iSequence] = value;
2539          double movement = 0.0;
2540
2541          if (value > tolerance) {
2542            // to lower bound (if swap)
2543            which[numberInfeasibilities++] = iSequence;
2544            movement = lower[iSequence] - upper[iSequence];
2545#ifndef NDEBUG
2546            if (fabs(movement) >= 1.0e30)
2547              resetFakeBounds(-1000 - iSequence);
2548#endif
2549#ifdef CLP_DEBUG
2550            if ((handler_->logLevel() & 32))
2551              printf("%d %d, new dj %g, alpha %g, movement %g\n",
2552                1, iSequence, value, alphaI, movement);
2553#endif
2554            changeObj += movement * cost[iSequence];
2555            matrix_->add(this, outputArray, iSequence, movement);
2556          }
2557        } else if (status == isFree) {
2558          double value = reducedCost[iSequence] - theta * alphaI;
2559          reducedCost[iSequence] = value;
2560        }
2561      }
2562    }
2563  } else {
2564    double *COIN_RESTRICT solution = solutionRegion(0);
2565    double *COIN_RESTRICT reducedCost = djRegion(0);
2566    double *COIN_RESTRICT lower = lowerRegion(0);
2567    double *COIN_RESTRICT upper = upperRegion(0);
2568    const double *COIN_RESTRICT cost = costRegion(0);
2569    int *COIN_RESTRICT which;
2570    which = rowArray->getIndices();
2571    int iSequence;
2572    for (iSequence = 0; iSequence < numberRows_; iSequence++) {
2573      double value = reducedCost[iSequence];
2574
2575      Status status = getStatus(iSequence + numberColumns_);
2576      // more likely to be at upper bound ?
2577      if (status == atUpperBound) {
2578        double movement = 0.0;
2579        //#define NO_SWAP7
2580        if (value > tolerance) {
2581          // to lower bound (if swap)
2582          // put back alpha
2583          which[numberInfeasibilities++] = iSequence;
2584          movement = lower[iSequence] - upper[iSequence];
2585#define TRY_SET_FAKE
2586#ifdef TRY_SET_FAKE
2587          if (fabs(movement) > dualBound_) {
2588            FakeBound bound = getFakeBound(iSequence + numberColumns_);
2589            if (bound == ClpSimplexDual::noFake) {
2590              setFakeBound(iSequence + numberColumns_,
2591                ClpSimplexDual::lowerFake);
2592              lower[iSequence] = upper[iSequence] - dualBound_;
2593              assert(fabs(lower[iSequence]) < 1.0e30);
2594              movement = lower[iSequence] - upper[iSequence];
2595              numberFake_++;
2596#ifndef NDEBUG
2597            } else {
2598              if (fabs(movement) >= 1.0e30)
2599                resetFakeBounds(-1000 - iSequence);
2600#endif
2601            }
2602          }
2603#endif
2604          changeObj += movement * cost[iSequence];
2605          outputArray->quickAdd(iSequence, -movement);
2606#ifndef NO_SWAP7
2607        } else if (value > -tolerance) {
2608          // at correct bound but may swap
2609          FakeBound bound = getFakeBound(iSequence + numberColumns_);
2610          if (bound == ClpSimplexDual::upperFake) {
2611            movement = lower[iSequence] - upper[iSequence];
2612#ifndef NDEBUG
2613            if (fabs(movement) >= 1.0e30)
2614              resetFakeBounds(-1000 - iSequence);
2615#endif
2616            setStatus(iSequence + numberColumns_, atLowerBound);
2617            solution[iSequence] = lower[iSequence];
2618            changeObj += movement * cost[iSequence];
2619            //numberFake_--;
2620            //setFakeBound(iSequence+numberColumns_,noFake);
2621          }
2622#endif
2623        }
2624      } else if (status == atLowerBound) {
2625        double movement = 0.0;
2626
2627        if (value < -tolerance) {
2628          // to upper bound
2629          // put back alpha
2630          which[numberInfeasibilities++] = iSequence;
2631          movement = upper[iSequence] - lower[iSequence];
2632#ifdef TRY_SET_FAKE
2633          if (fabs(movement) > dualBound_) {
2634            FakeBound bound = getFakeBound(iSequence + numberColumns_);
2635            if (bound == ClpSimplexDual::noFake) {
2636              setFakeBound(iSequence + numberColumns_,
2637                ClpSimplexDual::upperFake);
2638              upper[iSequence] = lower[iSequence] + dualBound_;
2639              assert(fabs(upper[iSequence]) < 1.0e30);
2640              movement = upper[iSequence] - lower[iSequence];
2641              numberFake_++;
2642#ifndef NDEBUG
2643            } else {
2644              if (fabs(movement) >= 1.0e30)
2645                resetFakeBounds(-1000 - iSequence);
2646#endif
2647            }
2648          }
2649#endif
2650          changeObj += movement * cost[iSequence];
2651          outputArray->quickAdd(iSequence, -movement);
2652#ifndef NO_SWAP7
2653        } else if (value < tolerance) {
2654          // at correct bound but may swap
2655          FakeBound bound = getFakeBound(iSequence + numberColumns_);
2656          if (bound == ClpSimplexDual::lowerFake) {
2657            movement = upper[iSequence] - lower[iSequence];
2658#ifndef NDEBUG
2659            if (fabs(movement) >= 1.0e30)
2660              resetFakeBounds(-1000 - iSequence);
2661#endif
2662            setStatus(iSequence + numberColumns_, atUpperBound);
2663            solution[iSequence] = upper[iSequence];
2664            changeObj += movement * cost[iSequence];
2665            //numberFake_--;
2666            //setFakeBound(iSequence+numberColumns_,noFake);
2667          }
2668#endif
2669        }
2670      }
2671    }
2672    // Do columns
2673    solution = solutionRegion(1);
2674    reducedCost = djRegion(1);
2675    lower = lowerRegion(1);
2676    upper = upperRegion(1);
2677    cost = costRegion(1);
2678    // set number of infeasibilities in row array
2679    numberRowInfeasibilities = numberInfeasibilities;
2680    rowArray->setNumElements(numberInfeasibilities);
2681    numberInfeasibilities = 0;
2682    which = columnArray->getIndices();
2683    for (iSequence = 0; iSequence < numberColumns_; iSequence++) {
2684      double value = reducedCost[iSequence];
2685
2686      Status status = getStatus(iSequence);
2687      if (status == atLowerBound) {
2688        double movement = 0.0;
2689
2690        if (value < -tolerance) {
2691          // to upper bound
2692          // put back alpha
2693          which[numberInfeasibilities++] = iSequence;
2694          movement = upper[iSequence] - lower[iSequence];
2695#ifdef TRY_SET_FAKE
2696          if (fabs(movement) > dualBound_) {
2697            FakeBound bound = getFakeBound(iSequence);
2698            if (bound == ClpSimplexDual::noFake) {
2699              setFakeBound(iSequence,
2700                ClpSimplexDual::upperFake);
2701              upper[iSequence] = lower[iSequence] + dualBound_;
2702              assert(fabs(upper[iSequence]) < 1.0e30);
2703              movement = upper[iSequence] - lower[iSequence];
2704              numberFake_++;
2705#ifndef NDEBUG
2706            } else {
2707              if (fabs(movement) >= 1.0e30)
2708                resetFakeBounds(-1000 - iSequence);
2709#endif
2710            }
2711          }
2712#endif
2713          changeObj += movement * cost[iSequence];
2714          matrix_->add(this, outputArray, iSequence, movement);
2715#ifndef NO_SWAP7
2716        } else if (value < tolerance) {
2717          // at correct bound but may swap
2718          FakeBound bound = getFakeBound(iSequence);
2719          if (bound == ClpSimplexDual::lowerFake) {
2720            movement = upper[iSequence] - lower[iSequence];
2721#ifndef NDEBUG
2722            if (fabs(movement) >= 1.0e30)
2723              resetFakeBounds(-1000 - iSequence);
2724#endif
2725            setStatus(iSequence, atUpperBound);
2726            solution[iSequence] = upper[iSequence];
2727            changeObj += movement * cost[iSequence];
2728            //numberFake_--;
2729            //setFakeBound(iSequence,noFake);
2730          }
2731#endif
2732        }
2733      } else if (status == atUpperBound) {
2734        double movement = 0.0;
2735
2736        if (value > tolerance) {
2737          // to lower bound (if swap)
2738          // put back alpha
2739          which[numberInfeasibilities++] = iSequence;
2740          movement = lower[iSequence] - upper[iSequence];
2741#ifdef TRY_SET_FAKE
2742          if (fabs(movement) > dualBound_) {
2743            FakeBound bound = getFakeBound(iSequence);
2744            if (bound == ClpSimplexDual::noFake) {
2745              setFakeBound(iSequence,
2746                ClpSimplexDual::lowerFake);
2747              lower[iSequence] = upper[iSequence] - dualBound_;
2748              assert(fabs(lower[iSequence]) < 1.0e30);
2749              movement = lower[iSequence] - upper[iSequence];
2750              numberFake_++;
2751#ifndef NDEBUG
2752            } else {
2753              if (fabs(movement) >= 1.0e30)
2754                resetFakeBounds(-1000 - iSequence);
2755#endif
2756            }
2757          }
2758#endif
2759          changeObj += movement * cost[iSequence];
2760          matrix_->add(this, outputArray, iSequence, movement);
2761#ifndef NO_SWAP7
2762        } else if (value > -tolerance) {
2763          // at correct bound but may swap
2764          FakeBound bound = getFakeBound(iSequence);
2765          if (bound == ClpSimplexDual::upperFake) {
2766            movement = lower[iSequence] - upper[iSequence];
2767#ifndef NDEBUG
2768            if (fabs(movement) >= 1.0e30)
2769              resetFakeBounds(-1000 - iSequence);
2770#endif
2771            setStatus(iSequence, atLowerBound);
2772            solution[iSequence] = lower[iSequence];
2773            changeObj += movement * cost[iSequence];
2774            //numberFake_--;
2775            //setFakeBound(iSequence,noFake);
2776          }
2777#endif
2778        }
2779      }
2780    }
2781  }
2782
2783#ifdef CLP_DEBUG
2784  if (fullRecompute && numberFake_ && (handler_->logLevel() & 16) != 0)
2785    printf("%d fake after full update\n", numberFake_);
2786#endif
2787  // set number of infeasibilities
2788  columnArray->setNumElements(numberInfeasibilities);
2789  numberInfeasibilities += numberRowInfeasibilities;
2790  if (fullRecompute) {
2791    // do actual flips
2792    flipBounds(rowArray, columnArray);
2793  }
2794  objectiveChange += changeObj;
2795  return numberInfeasibilities;
2796}
2797void ClpSimplexDual::updateDualsInValuesPass(CoinIndexedVector *rowArray,
2798  CoinIndexedVector *columnArray,
2799  double theta)
2800{
2801
2802  // use a tighter tolerance except for all being okay
2803  double tolerance = dualTolerance_;
2804
2805  // Coding is very similar but we can save a bit by splitting
2806  // Do rows
2807  {
2808    int i;
2809    double *reducedCost = djRegion(0);
2810    double *work;
2811    int number;
2812    int *which;
2813    work = rowArray->denseVector();
2814    number = rowArray->getNumElements();
2815    which = rowArray->getIndices();
2816    for (i = 0; i < number; i++) {
2817      int iSequence = which[i];
2818      double alphaI = work[i];
2819      double value = reducedCost[iSequence] - theta * alphaI;
2820      work[i] = 0.0;
2821      reducedCost[iSequence] = value;
2822
2823      Status status = getStatus(iSequence + numberColumns_);
2824      // more likely to be at upper bound ?
2825      if (status == atUpperBound) {
2826
2827        if (value > tolerance)
2828          reducedCost[iSequence] = 0.0;
2829      } else if (status == atLowerBound) {
2830
2831        if (value < -tolerance) {
2832          reducedCost[iSequence] = 0.0;
2833        }
2834      }
2835    }
2836  }
2837  rowArray->setNumElements(0);
2838
2839  // Do columns
2840  {
2841    int i;
2842    double *reducedCost = djRegion(1);
2843    double *work;
2844    int number;
2845    int *which;
2846    work = columnArray->denseVector();
2847    number = columnArray->getNumElements();
2848    which = columnArray->getIndices();
2849
2850    for (i = 0; i < number; i++) {
2851      int iSequence = which[i];
2852      double alphaI = work[i];
2853      double value = reducedCost[iSequence] - theta * alphaI;
2854      work[i] = 0.0;
2855      reducedCost[iSequence] = value;
2856
2857      Status status = getStatus(iSequence);
2858      if (status == atLowerBound) {
2859        if (value < -tolerance)
2860          reducedCost[iSequence] = 0.0;
2861      } else if (status == atUpperBound) {
2862        if (value > tolerance)
2863          reducedCost[iSequence] = 0.0;
2864      }
2865    }
2866  }
2867  columnArray->setNumElements(0);
2868}
2869/*
2870   Chooses dual pivot row
2871   Would be faster with separate region to scan
2872   and will have this (with square of infeasibility) when steepest
2873   For easy problems we can just choose one of the first rows we look at
2874*/
2875void ClpSimplexDual::dualRow(int alreadyChosen)
2876{
2877  // get pivot row using whichever method it is
2878  int chosenRow = -1;
2879#ifdef FORCE_FOLLOW
2880  bool forceThis = false;
2881  if (!fpFollow && strlen(forceFile)) {
2882    fpFollow = fopen(forceFile, "r");
2883    assert(fpFollow);
2884  }
2885  if (fpFollow) {
2886    if (numberIterations_ <= force_iteration) {
2887      // read to next Clp0102
2888      char temp[300];
2889      while (fgets(temp, 250, fpFollow)) {
2890        if (strncmp(temp, "Clp0102", 7))
2891          continue;
2892        char cin, cout;
2893        sscanf(temp + 9, "%d%*f%*s%*c%c%d%*s%*c%c%d",
2894          &force_iteration, &cin, &force_in, &cout, &force_out);
2895        if (cin == 'R')
2896          force_in += numberColumns_;
2897        if (cout == 'R')
2898          force_out += numberColumns_;
2899        forceThis = true;
2900        assert(numberIterations_ == force_iteration - 1);
2901        printf("Iteration %d will force %d out and %d in\n",
2902          force_iteration, force_out, force_in);
2903        alreadyChosen = force_out;
2904        break;
2905      }
2906    } else {
2907      // use old
2908      forceThis = true;
2909    }
2910    if (!forceThis) {
2911      fclose(fpFollow);
2912      fpFollow = NULL;
2913      forceFile = "";
2914    }
2915  }
2916#endif
2917  //double freeAlpha = 0.0;
2918  if (alreadyChosen < 0) {
2919    // first see if any free variables and put them in basis
2920    int nextFree = nextSuperBasic();
2921    //nextFree=-1; //off
2922    if (nextFree >= 0) {
2923      // unpack vector and find a good pivot
2924      unpack(rowArray_[1], nextFree);
2925      factorization_->updateColumn(rowArray_[2], rowArray_[1]);
2926
2927      double *work = rowArray_[1]->denseVector();
2928      int number = rowArray_[1]->getNumElements();
2929      int *which = rowArray_[1]->getIndices();
2930      double bestFeasibleAlpha = 0.0;
2931      int bestFeasibleRow = -1;
2932      double bestInfeasibleAlpha = 0.0;
2933      int bestInfeasibleRow = -1;
2934      int i;
2935
2936      for (i = 0; i < number; i++) {
2937        int iRow = which[i];
2938        double alpha = fabs(work[iRow]);
2939        if (alpha > 1.0e-3) {
2940          int iSequence = pivotVariable_[iRow];
2941          double value = solution_[iSequence];
2942          double lower = lower_[iSequence];
2943          double upper = upper_[iSequence];
2944          double infeasibility = 0.0;
2945          if (value > upper)
2946            infeasibility = value - upper;
2947          else if (value < lower)
2948            infeasibility = lower - value;
2949          if (infeasibility * alpha > bestInfeasibleAlpha && alpha > 1.0e-1) {
2950            if (!flagged(iSequence)) {
2951              bestInfeasibleAlpha = infeasibility * alpha;
2952              bestInfeasibleRow = iRow;
2953            }
2954          }
2955          if (alpha > bestFeasibleAlpha && (lower > -1.0e20 || upper < 1.0e20)) {
2956            bestFeasibleAlpha = alpha;
2957            bestFeasibleRow = iRow;
2958          }
2959        }
2960      }
2961      if (bestInfeasibleRow >= 0)
2962        chosenRow = bestInfeasibleRow;
2963      else if (bestFeasibleAlpha > 1.0e-2)
2964        chosenRow = bestFeasibleRow;
2965      if (chosenRow >= 0) {
2966        pivotRow_ = chosenRow;
2967        //freeAlpha = work[chosenRow];
2968      }
2969      rowArray_[1]->clear();
2970    }
2971  } else {
2972    // in values pass
2973    chosenRow = alreadyChosen;
2974#ifdef FORCE_FOLLOW
2975    if (forceThis) {
2976      alreadyChosen = -1;
2977      chosenRow = -1;
2978      for (int i = 0; i < numberRows_; i++) {
2979        if (pivotVariable_[i] == force_out) {
2980          chosenRow = i;
2981          break;
2982        }
2983      }
2984      assert(chosenRow >= 0);
2985    }
2986#endif
2987    pivotRow_ = chosenRow;
2988  }
2989  if (chosenRow < 0)
2990    pivotRow_ = dualRowPivot_->pivotRow();
2991
2992  if (pivotRow_ >= 0) {
2993    sequenceOut_ = pivotVariable_[pivotRow_];
2994    valueOut_ = solution_[sequenceOut_];
2995    lowerOut_ = lower_[sequenceOut_];
2996    upperOut_ = upper_[sequenceOut_];
2997    if (alreadyChosen < 0) {
2998      // if we have problems we could try other way and hope we get a
2999      // zero pivot?
3000      if (valueOut_ > upperOut_) {
3001        directionOut_ = -1;
3002        dualOut_ = valueOut_ - upperOut_;
3003      } else if (valueOut_ < lowerOut_) {
3004        directionOut_ = 1;
3005        dualOut_ = lowerOut_ - valueOut_;
3006      } else {
3007#if 1
3008        // odd (could be free) - it's feasible - go to nearest
3009        if (valueOut_ - lowerOut_ < upperOut_ - valueOut_) {
3010          directionOut_ = 1;
3011          dualOut_ = lowerOut_ - valueOut_;
3012        } else {
3013          directionOut_ = -1;
3014          dualOut_ = valueOut_ - upperOut_;
3015        }
3016#else
3017        // odd (could be free) - it's feasible - improve obj
3018        printf("direction from alpha of %g is %d\n",
3019          freeAlpha, freeAlpha > 0.0 ? 1 : -1);
3020        if (valueOut_ - lowerOut_ > 1.0e20)
3021          freeAlpha = 1.0;
3022        else if (upperOut_ - valueOut_ > 1.0e20)
3023          freeAlpha = -1.0;
3024        //if (valueOut_-lowerOut_<upperOut_-valueOut_) {
3025        if (freeAlpha < 0.0) {
3026          directionOut_ = 1;
3027          dualOut_ = lowerOut_ - valueOut_;
3028        } else {
3029          directionOut_ = -1;
3030          dualOut_ = valueOut_ - upperOut_;
3031        }
3032        printf("direction taken %d - bounds %g %g %g\n",
3033          directionOut_, lowerOut_, valueOut_, upperOut_);
3034#endif
3035      }
3036#ifdef CLP_DEBUG
3037      assert(dualOut_ >= 0.0);
3038#endif
3039    } else {
3040      // in values pass so just use sign of dj
3041      // We don't want to go through any barriers so set dualOut low
3042      // free variables will never be here
3043      dualOut_ = 1.0e-6;
3044      if (dj_[sequenceOut_] > 0.0) {
3045        // this will give a -1 in pivot row (as slacks are -1.0)
3046        directionOut_ = 1;
3047      } else {
3048        directionOut_ = -1;
3049      }
3050    }
3051  }
3052  return;
3053}
3054// Checks if any fake bounds active - if so returns number and modifies
3055// dualBound_ and everything.
3056// Free variables will be left as free
3057// Returns number of bounds changed if >=0
3058// Returns -1 if not initialize and no effect
3059// Fills in changeVector which can be used to see if unbounded
3060// and cost of change vector
3061int ClpSimplexDual::changeBounds(int initialize,
3062  CoinIndexedVector *outputArray,
3063  double &changeCost)
3064{
3065  numberFake_ = 0;
3066  if (!initialize) {
3067    int numberInfeasibilities;
3068    double newBound;
3069    newBound = 5.0 * dualBound_;
3070    numberInfeasibilities = 0;
3071    changeCost = 0.0;
3072    // put back original bounds and then check
3073    createRim1(false);
3074    int iSequence;
3075    // bounds will get bigger - just look at ones at bounds
3076    for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
3077      double lowerValue = lower_[iSequence];
3078      double upperValue = upper_[iSequence];
3079      double value = solution_[iSequence];
3080      setFakeBound(iSequence, ClpSimplexDual::noFake);
3081      switch (getStatus(iSequence)) {
3082
3083      case basic:
3084      case ClpSimplex::isFixed:
3085        break;
3086      case isFree:
3087      case superBasic:
3088        break;
3089      case atUpperBound:
3090        if (fabs(value - upperValue) > primalTolerance_) {
3091          if (fabs(dj_[iSequence]) > 1.0e-9) {
3092            numberInfeasibilities++;
3093          } else {
3094            setStatus(iSequence, superBasic);
3095            moreSpecialOptions_ &= ~8;
3096          }
3097        }
3098        break;
3099      case atLowerBound:
3100        if (fabs(value - lowerValue) > primalTolerance_) {
3101          if (fabs(dj_[iSequence]) > 1.0e-9) {
3102            numberInfeasibilities++;
3103          } else {
3104            setStatus(iSequence, superBasic);
3105            moreSpecialOptions_ &= ~8;
3106          }
3107        }
3108        break;
3109      }
3110    }
3111    // If dual infeasible then carry on
3112    if (numberInfeasibilities) {
3113      handler_->message(CLP_DUAL_CHECKB, messages_)
3114        << newBound
3115        << CoinMessageEol;
3116      int iSequence;
3117      for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
3118        double lowerValue = lower_[iSequence];
3119        double upperValue = upper_[iSequence];
3120        double newLowerValue;
3121        double newUpperValue;
3122        Status status = getStatus(iSequence);
3123        if (status == atUpperBound || status == atLowerBound) {
3124          double value = solution_[iSequence];
3125          if (value - lowerValue <= upperValue - value) {
3126            newLowerValue = CoinMax(lowerValue, value - 0.666667 * newBound);
3127            newUpperValue = CoinMin(upperValue, newLowerValue + newBound);
3128          } else {
3129            newUpperValue = CoinMin(upperValue, value + 0.666667 * newBound);
3130            newLowerValue = CoinMax(lowerValue, newUpperValue - newBound);
3131          }
3132          if (newLowerValue > lowerValue) {
3133            if (newUpperValue < upperValue) {
3134              setFakeBound(iSequence, ClpSimplexDual::bothFake);
3135              // redo
3136              if (status == atLowerBound) {
3137                newLowerValue = value;
3138                newUpperValue = CoinMin(upperValue, newLowerValue + newBound);
3139              } else {
3140                newUpperValue = value;
3141                newLowerValue = CoinMax(lowerValue, newUpperValue - newBound);
3142              }
3143              numberFake_++;
3144            } else {
3145              setFakeBound(iSequence, ClpSimplexDual::lowerFake);
3146              numberFake_++;
3147            }
3148          } else {
3149            if (newUpperValue < upperValue) {
3150              setFakeBound(iSequence, ClpSimplexDual::upperFake);
3151              numberFake_++;
3152            }
3153          }
3154          lower_[iSequence] = newLowerValue;
3155          upper_[iSequence] = newUpperValue;
3156          if (status == atUpperBound)
3157            solution_[iSequence] = newUpperValue;
3158          else
3159            solution_[iSequence] = newLowerValue;
3160          double movement = solution_[iSequence] - value;
3161          if (movement && outputArray) {
3162            if (iSequence >= numberColumns_) {
3163              outputArray->quickAdd(iSequence, -movement);
3164              changeCost += movement * cost_[iSequence];
3165            } else {
3166              matrix_->add(this, outputArray, iSequence, movement);
3167              changeCost += movement * cost_[iSequence];
3168            }
3169          }
3170        }
3171      }
3172      dualBound_ = newBound;
3173    } else {
3174      numberInfeasibilities = -1;
3175    }
3176    return numberInfeasibilities;
3177  } else if (initialize == 1 || initialize == 3) {
3178    int iSequence;
3179    if (initialize == 3) {
3180      if (columnScale_) {
3181        for (iSequence = 0; iSequence < numberColumns_; iSequence++) {
3182          if (getFakeBound(iSequence) != ClpSimplexDual::noFake) {
3183            double multiplier = rhsScale_ * inverseColumnScale_[iSequence];
3184            // lower
3185            double value = columnLower_[iSequence];
3186            if (value > -1.0e30) {
3187              value *= multiplier;
3188            }
3189            lower_[iSequence] = value;
3190            // upper
3191            value = columnUpper_[iSequence];
3192            if (value < 1.0e30) {
3193              value *= multiplier;
3194            }
3195            upper_[iSequence] = value;
3196            setFakeBound(iSequence, ClpSimplexDual::noFake);
3197          }
3198        }
3199        for (iSequence = 0; iSequence < numberRows_; iSequence++) {
3200          // lower
3201          double multiplier = rhsScale_ * rowScale_[iSequence];
3202          double value = rowLower_[iSequence];
3203          if (value > -1.0e30) {
3204            value *= multiplier;
3205          }
3206          lower_[iSequence + numberColumns_] = value;
3207          // upper
3208          value = rowUpper_[iSequence];
3209          if (value < 1.0e30) {
3210            value *= multiplier;
3211          }
3212          upper_[iSequence + numberColumns_] = value;
3213          setFakeBound(iSequence + numberColumns_, ClpSimplexDual::noFake);
3214        }
3215      } else {
3216        for (iSequence = 0; iSequence < numberColumns_; iSequence++) {
3217          if (getFakeBound(iSequence) != ClpSimplexDual::noFake) {
3218            lower_[iSequence] = columnLower_[iSequence];
3219            upper_[iSequence] = columnUpper_[iSequence];
3220            setFakeBound(iSequence, ClpSimplexDual::noFake);
3221          }
3222        }
3223        for (iSequence = 0; iSequence < numberRows_; iSequence++) {
3224          if (getFakeBound(iSequence + numberColumns_) != ClpSimplexDual::noFake) {
3225            lower_[iSequence + numberColumns_] = rowLower_[iSequence];
3226            upper_[iSequence + numberColumns_] = rowUpper_[iSequence];
3227            setFakeBound(iSequence + numberColumns_, ClpSimplexDual::noFake);
3228          }
3229        }
3230      }
3231    }
3232    double testBound = 0.999999 * dualBound_;
3233    for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
3234      Status status = getStatus(iSequence);
3235      if (status == atUpperBound || status == atLowerBound) {
3236        double lowerValue = lower_[iSequence];
3237        double upperValue = upper_[iSequence];
3238        double value = solution_[iSequence];
3239        if (lowerValue > -largeValue_ || upperValue < largeValue_) {
3240          if (true || lowerValue - value > -0.5 * dualBound_ || upperValue - value < 0.5 * dualBound_) {
3241            if (fabs(lowerValue - value) <= fabs(upperValue - value)) {
3242              if (upperValue > lowerValue + testBound) {
3243                if (getFakeBound(iSequence) == ClpSimplexDual::noFake)
3244                  numberFake_++;
3245                upper_[iSequence] = lowerValue + dualBound_;
3246                setFakeBound(iSequence, ClpSimplexDual::upperFake);
3247              }
3248            } else {
3249              if (lowerValue < upperValue - testBound) {
3250                if (getFakeBound(iSequence) == ClpSimplexDual::noFake)
3251                  numberFake_++;
3252                lower_[iSequence] = upperValue - dualBound_;
3253                setFakeBound(iSequence, ClpSimplexDual::lowerFake);
3254              }
3255            }
3256          } else {
3257            if (getFakeBound(iSequence) == ClpSimplexDual::noFake)
3258              numberFake_++;
3259            lower_[iSequence] = -0.5 * dualBound_;
3260            upper_[iSequence] = 0.5 * dualBound_;
3261            setFakeBound(iSequence, ClpSimplexDual::bothFake);
3262            abort();
3263          }
3264          if (status == atUpperBound)
3265            solution_[iSequence] = upper_[iSequence];
3266          else
3267            solution_[iSequence] = lower_[iSequence];
3268        } else {
3269          // set non basic free variables to fake bounds
3270          // I don't think we should ever get here
3271          // yes we can if basis goes singular twice in succession!
3272          //CoinAssert(!("should not be here"));
3273          lower_[iSequence] = -0.5 * dualBound_;
3274          upper_[iSequence] = 0.5 * dualBound_;
3275          setFakeBound(iSequence, ClpSimplexDual::bothFake);
3276          numberFake_++;
3277          setStatus(iSequence, atUpperBound);
3278          solution_[iSequence] = 0.5 * dualBound_;
3279        }
3280      } else if (status == basic) {
3281        // make sure not at fake bound and bounds correct
3282        setFakeBound(iSequence, ClpSimplexDual::noFake);
3283        double gap = upper_[iSequence] - lower_[iSequence];
3284        if (gap > 0.5 * dualBound_ && gap < 2.0 * dualBound_) {
3285          if (iSequence < numberColumns_) {
3286            if (columnScale_) {
3287              double multiplier = rhsScale_ * inverseColumnScale_[iSequence];
3288              // lower
3289              double value = columnLower_[iSequence];
3290              if (value > -1.0e30) {
3291                value *= multiplier;
3292              }
3293              lower_[iSequence] = value;
3294              // upper
3295              value = columnUpper_[iSequence];
3296              if (value < 1.0e30) {
3297                value *= multiplier;
3298              }
3299              upper_[iSequence] = value;
3300            } else {
3301              lower_[iSequence] = columnLower_[iSequence];
3302              ;
3303              upper_[iSequence] = columnUpper_[iSequence];
3304              ;
3305            }
3306          } else {
3307            int iRow = iSequence - numberColumns_;
3308            if (rowScale_) {
3309              // lower
3310              double multiplier = rhsScale_ * rowScale_[iRow];
3311              double value = rowLower_[iRow];
3312              if (value > -1.0e30) {
3313                value *= multiplier;
3314              }
3315              lower_[iSequence] = value;
3316              // upper
3317              value = rowUpper_[iRow];
3318              if (value < 1.0e30) {
3319                value *= multiplier;
3320              }
3321              upper_[iSequence] = value;
3322            } else {
3323              lower_[iSequence] = rowLower_[iRow];
3324              ;
3325              upper_[iSequence] = rowUpper_[iRow];
3326              ;
3327            }
3328          }
3329        }
3330      }
3331    }
3332
3333    return 1;
3334  } else {
3335    // just reset changed ones
3336    if (columnScale_) {
3337      int iSequence;
3338      for (iSequence = 0; iSequence < numberColumns_; iSequence++) {
3339        FakeBound fakeStatus = getFakeBound(iSequence);
3340        if (fakeStatus != noFake) {
3341          if ((static_cast< int >(fakeStatus) & 1) != 0) {
3342            // lower
3343            double value = columnLower_[iSequence];
3344            if (value > -1.0e30) {
3345              double multiplier = rhsScale_ * inverseColumnScale_[iSequence];
3346              value *= multiplier;
3347            }
3348            columnLowerWork_[iSequence] = value;
3349          }
3350          if ((static_cast< int >(fakeStatus) & 2) != 0) {
3351            // upper
3352            double value = columnUpper_[iSequence];
3353            if (value < 1.0e30) {
3354              double multiplier = rhsScale_ * inverseColumnScale_[iSequence];
3355              value *= multiplier;
3356            }
3357            columnUpperWork_[iSequence] = value;
3358          }
3359        }
3360      }
3361      for (iSequence = 0; iSequence < numberRows_; iSequence++) {
3362        FakeBound fakeStatus = getFakeBound(iSequence + numberColumns_);
3363        if (fakeStatus != noFake) {
3364          if ((static_cast< int >(fakeStatus) & 1) != 0) {
3365            // lower
3366            double value = rowLower_[iSequence];
3367            if (value > -1.0e30) {
3368              double multiplier = rhsScale_ * rowScale_[iSequence];
3369              value *= multiplier;
3370            }
3371            rowLowerWork_[iSequence] = value;
3372          }
3373          if ((static_cast< int >(fakeStatus) & 2) != 0) {
3374            // upper
3375            double value = rowUpper_[iSequence];
3376            if (value < 1.0e30) {
3377              double multiplier = rhsScale_ * rowScale_[iSequence];
3378              value *= multiplier;
3379            }
3380            rowUpperWork_[iSequence] = value;
3381          }
3382        }
3383      }
3384    } else {
3385      int iSequence;
3386      for (iSequence = 0; iSequence < numberColumns_; iSequence++) {
3387        FakeBound fakeStatus = getFakeBound(iSequence);
3388        if ((static_cast< int >(fakeStatus) & 1) != 0) {
3389          // lower
3390          columnLowerWork_[iSequence] = columnLower_[iSequence];
3391        }
3392        if ((static_cast< int >(fakeStatus) & 2) != 0) {
3393          // upper
3394          columnUpperWork_[iSequence] = columnUpper_[iSequence];
3395        }
3396      }
3397      for (iSequence = 0; iSequence < numberRows_; iSequence++) {
3398        FakeBound fakeStatus = getFakeBound(iSequence + numberColumns_);
3399        if ((static_cast< int >(fakeStatus) & 1) != 0) {
3400          // lower
3401          rowLowerWork_[iSequence] = rowLower_[iSequence];
3402        }
3403        if ((static_cast< int >(fakeStatus) & 2) != 0) {
3404          // upper
3405          rowUpperWork_[iSequence] = rowUpper_[iSequence];
3406        }
3407      }
3408    }
3409    return 0;
3410  }
3411}
3412// Just checks if any fake bounds active - if so returns number
3413int ClpSimplexDual::checkFakeBounds() const
3414{
3415  int numberActive = 0;
3416  for (int iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
3417    switch (getStatus(iSequence)) {
3418
3419    case basic:
3420    case ClpSimplex::isFixed:
3421      break;
3422    case isFree:
3423    case superBasic:
3424      break;
3425    case atUpperBound:
3426      if ((getFakeBound(iSequence) & 2) != 0)
3427        numberActive++;
3428      break;
3429    case atLowerBound:
3430      if ((getFakeBound(iSequence) & 1) != 0)
3431        numberActive++;
3432      break;
3433    }
3434  }
3435  return numberActive;
3436}
3437#if ABOCA_LITE
3438/* Meat of transposeTimes by column when not scaled and skipping
3439   and doing part of dualColumn */
3440static void
3441dualColumn00(clpTempInfo &info)
3442{
3443  const int *COIN_RESTRICT which = info.which;
3444  const double *COIN_RESTRICT work = info.work;
3445  int *COIN_RESTRICT index = info.index;
3446  double *COIN_RESTRICT spare = info.spare;
3447  const unsigned char *COIN_RESTRICT status = info.status;
3448  const double *COIN_RESTRICT reducedCost = info.reducedCost;
3449  double upperTheta = info.upperTheta;
3450  double acceptablePivot = info.acceptablePivot;
3451  double dualTolerance = info.tolerance;
3452  int numberToDo = info.numberToDo;
3453  double tentativeTheta = 1.0e15;
3454  int numberRemaining = 0;
3455  double multiplier[] = { -1.0, 1.0 };
3456  double dualT = -dualTolerance;
3457  for (int i = 0; i < numberToDo; i++) {
3458    int iSequence = which[i];
3459    int wanted = (status[iSequence] & 3) - 1;
3460    if (wanted) {
3461      double mult = multiplier[wanted - 1];
3462      double alpha = work[i] * mult;
3463      if (alpha > 0.0) {
3464        double oldValue = reducedCost[iSequence] * mult;
3465        double value = oldValue - tentativeTheta * alpha;
3466        if (value < dualT) {
3467          value = oldValue - upperTheta * alpha;
3468          if (value < dualT && alpha >= acceptablePivot) {
3469            upperTheta = (oldValue - dualT) / alpha;
3470          }
3471          // add to list
3472          spare[numberRemaining] = alpha * mult;
3473          index[numberRemaining++] = iSequence;
3474        }
3475      }
3476    }
3477  }
3478  info.numberRemaining = numberRemaining;
3479  info.upperTheta = upperTheta;
3480}
3481static void
3482dualColumn000(int numberThreads, clpTempInfo *info)
3483{
3484  for (int i = 0; i < numberThreads; i++) {
3485    cilk_spawn dualColumn00(info[i]);
3486  }
3487  cilk_sync;
3488}
3489void moveAndZero(clpTempInfo *info, int type, void *extra)
3490{
3491  int numberThreads = abcState();
3492  switch (type) {
3493  case 1: {
3494    int numberRemaining = info[0].numberRemaining;
3495    int *COIN_RESTRICT index = info[0].index + numberRemaining;
3496    double *COIN_RESTRICT spare = info[0].spare + numberRemaining;
3497    for (int i = 1; i < numberThreads; i++) {
3498      int number = info[i].numberRemaining;
3499      memmove(index, info[i].index, number * sizeof(int));
3500      index += number;
3501      double *COIN_RESTRICT from = info[i].spare;
3502      assert(from >= spare);
3503      memmove(spare, from, number * sizeof(double));
3504      spare += number;
3505    }
3506    // now zero out
3507    int i;
3508    for (i = 1; i < numberThreads; i++) {
3509      double *spareBit = info[i].spare + info[i].numberRemaining;
3510      if (spareBit > spare) {
3511        memset(spare, 0, (spareBit - spare) * sizeof(double));
3512        break;
3513      }
3514    }
3515    i++; // just zero
3516    for (; i < numberThreads; i++) {
3517      int number = info[i].numberRemaining;
3518      memset(info[i].spare, 0, number * sizeof(double));
3519    }
3520  } break;
3521  case 2: {
3522    int numberAdded = info[0].numberAdded;
3523    int *COIN_RESTRICT index = info[0].which + numberAdded;
3524    double *COIN_RESTRICT spare = info[0].infeas + numberAdded;
3525    for (int i = 1; i < numberThreads; i++) {
3526      int number = info[i].numberAdded;
3527      memmove(index, info[i].which, number * sizeof(int));
3528      index += number;
3529      double *COIN_RESTRICT from = info[i].infeas;
3530      assert(from >= spare);
3531      memmove(spare, from, number * sizeof(double));
3532      spare += number;
3533    }
3534    // now zero out
3535    int i;
3536    for (i = 1; i < numberThreads; i++) {
3537      double *spareBit = info[i].infeas + info[i].numberAdded;
3538      if (spareBit > spare) {
3539        memset(spare, 0, (spareBit - spare) * sizeof(double));
3540        break;
3541      }
3542    }
3543    i++; // just zero
3544    for (; i < numberThreads; i++) {
3545      int number = info[i].numberAdded;
3546      memset(info[i].infeas, 0, number * sizeof(double));
3547    }
3548  } break;
3549  default:
3550    abort();
3551    break;
3552  }
3553}
3554#endif
3555#ifdef _MSC_VER
3556#include <intrin.h>
3557#else
3558#include <immintrin.h>
3559//#include <fmaintrin.h>
3560#endif
3561int ClpSimplexDual::dualColumn0(const CoinIndexedVector *rowArray,
3562  const CoinIndexedVector *columnArray,
3563  CoinIndexedVector *spareArray,
3564  double acceptablePivot,
3565  double &upperReturn, double &badFree)
3566{
3567  // do first pass to get possibles
3568  double *spare = spareArray->denseVector();
3569  int *index = spareArray->getIndices();
3570  const double *work;
3571  int number;
3572  const int *which;
3573  const double *reducedCost;
3574  // We can also see if infeasible or pivoting on free
3575  double tentativeTheta = 1.0e15;
3576  double upperTheta = 1.0e31;
3577  double freePivot = acceptablePivot;
3578  int numberRemaining = 0;
3579  int i;
3580  badFree = 0.0;
3581  if ((moreSpecialOptions_ & 8) != 0) {
3582    // No free or super basic
3583    // bestPossible will re recomputed if necessary
3584#ifndef COIN_AVX2
3585    double multiplier[] = { -1.0, 1.0 };
3586#else
3587    double multiplier[4] = { 0.0, 0.0, -1.0, 1.0 };
3588#endif
3589    double dualT = -dualTolerance_;
3590#if ABOCA_LITE == 0
3591    int nSections = 2;
3592#else
3593    int numberThreads = abcState();
3594    int nSections = numberThreads ? 1 : 2;
3595#endif
3596    for (int iSection = 0; iSection < nSections; iSection++) {
3597
3598      int addSequence;
3599      unsigned char *statusArray;
3600      if (!iSection) {
3601        work = rowArray->denseVector();
3602        number = rowArray->getNumElements();
3603        which = rowArray->getIndices();
3604        reducedCost = rowReducedCost_;
3605        addSequence = numberColumns_;
3606        statusArray = status_ + numberColumns_;
3607      } else {
3608        work = columnArray->denseVector();
3609        number = columnArray->getNumElements();
3610        which = columnArray->getIndices();
3611        reducedCost = reducedCostWork_;
3612        addSequence = 0;
3613        statusArray = status_;
3614      }
3615#ifndef COIN_AVX2
3616      for (i = 0; i < number; i++) {
3617        int iSequence = which[i];
3618        double alpha;
3619        double oldValue;
3620        double value;
3621
3622        assert(getStatus(iSequence + addSequence) != isFree
3623          && getStatus(iSequence + addSequence) != superBasic);
3624        int iStatus = (statusArray[iSequence] & 3) - 1;
3625        if (iStatus) {
3626          double mult = multiplier[iStatus - 1];
3627          alpha = work[i] * mult;
3628          if (alpha > 0.0) {
3629            oldValue = reducedCost[iSequence] * mult;
3630            value = oldValue - tentativeTheta * alpha;
3631            if (value < dualT) {
3632              value = oldValue - upperTheta * alpha;
3633              if (value < dualT && alpha >= acceptablePivot) {
3634                upperTheta = (oldValue - dualT) / alpha;
3635                //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
3636              }
3637              // add to list
3638              spare[numberRemaining] = alpha * mult;
3639              index[numberRemaining++] = iSequence + addSequence;
3640            }
3641          }
3642        }
3643      }
3644      //
3645#else
3646      //#define COIN_AVX2 4 // temp
3647#if COIN_AVX2 == 1
3648#define COIN_AVX2_SHIFT 0
3649#elif COIN_AVX2 == 2
3650#define COIN_AVX2_SHIFT 1
3651#elif COIN_AVX2 == 4
3652#define COIN_AVX2_SHIFT 2
3653#elif COIN_AVX2 == 8
3654#define COIN_AVX2_SHIFT 3
3655#else
3656      error;
3657#endif
3658      //#define COIN_ALIGN 8*COIN_AVX2 // later
3659      //#define COIN_ALIGN_DOUBLE COIN_AVX2
3660#define CHECK_CHUNK 4
3661      // round up
3662      int *whichX = const_cast< int * >(which);
3663      double *workX = const_cast< double * >(work);
3664      int nBlocks = (number + CHECK_CHUNK - 1) / CHECK_CHUNK;
3665      int n = nBlocks * CHECK_CHUNK + 1;
3666      for (int i = number; i < n; i++) {
3667        workX[i] = 0.0;
3668        whichX[i] = 0; // alpha will be zero so not chosen
3669      }
3670      bool acceptableX[CHECK_CHUNK + 1];
3671      double oldValueX[CHECK_CHUNK + 1];
3672      double newValueX[CHECK_CHUNK + 1];
3673      double alphaX[CHECK_CHUNK + 1];
3674      newValueX[CHECK_CHUNK] = 0.0;
3675#define USE_USE_AVX
3676      //#define CHECK_H 1
3677#ifdef USE_USE_AVX
3678#define NEED_AVX
3679#elif CHECK_H
3680#define NEED_AVX
3681#endif
3682#ifdef NEED_AVX
3683      double mult2[CHECK_CHUNK] __attribute__((aligned(64)));
3684      CoinInt64 acceptableY[CHECK_CHUNK] __attribute__((aligned(64)));
3685      CoinInt64 goodDj[CHECK_CHUNK + 1] __attribute__((aligned(64)));
3686      double oldValueY[CHECK_CHUNK] __attribute__((aligned(64)));
3687      double alphaY[CHECK_CHUNK + 1] __attribute__((aligned(64)));
3688      memset(acceptableY, 0, sizeof(acceptableY));
3689      memset(goodDj, 0, sizeof(goodDj));
3690      memset(oldValueY, 0, sizeof(oldValueY));
3691      memset(alphaY, 0, sizeof(alphaY));
3692      __m256d tentative2 = _mm256_set1_pd(-tentativeTheta);
3693      __m256d dualT2 = _mm256_set1_pd(dualT);
3694      __m256d acceptable2 = _mm256_set1_pd(acceptablePivot);
3695#endif
3696      for (int iBlock = 0; iBlock < nBlocks; iBlock++) {
3697        bool store = false;
3698        double alpha = 0.0;
3699        double oldValue = 0.0;
3700        double newValue = 0.0;
3701        double trueAlpha = 0.0;
3702        int jSequence = 0;
3703#ifndef USE_USE_AVX
3704        for (int i = 0; i < CHECK_CHUNK + 1; i++) {
3705          int iSequence = which[i];
3706          int iStatus = (statusArray[iSequence] & 3);
3707          double mult = multiplier[iStatus];
3708          double newAlpha = work[i] * mult;
3709          double oldDj = reducedCost[iSequence] * mult;
3710          newValue = (oldDj - tentativeTheta * newAlpha) - dualT;
3711          acceptableX[i] = newAlpha >= acceptablePivot;
3712          oldValueX[i] = oldDj;
3713          newValueX[i] = newValue;
3714          alphaX[i] = newAlpha;
3715        }
3716#endif
3717#ifdef NEED_AVX
3718        __m128i columns = _mm_load_si128((const __m128i *)which);
3719        // what do we get - this must be wrong
3720        // probably only 1 and 2 - can we be clever
3721        // fix
3722        //__m128i status; // = _mm256_i32gather_ps(statusArray,columns,1);
3723        //status.m128i_i32[0]=statusArray[columns.m128i_i32[0]];
3724        for (int i = 0; i < CHECK_CHUNK; i++) {
3725          int iSequence = which[i];
3726          int iStatus = (statusArray[iSequence] & 3);
3727          mult2[i] = multiplier[iStatus];
3728        }
3729        //__m256d newAlpha2 = _mm256_i32gather_pd(multiplier,status,1); // mult here
3730        __m256d newAlpha2 = _mm256_load_pd(mult2);
3731        __m256d oldDj = _mm256_i32gather_pd(reducedCost, columns, 8);
3732        oldDj = _mm256_mul_pd(oldDj, newAlpha2); // remember newAlpha==mult
3733        _mm256_store_pd(oldValueY, oldDj); // redo later
3734        __m256d work2 = _mm256_load_pd(work);
3735        newAlpha2 = _mm256_mul_pd(newAlpha2, work2); // now really newAlpha
3736        //__m256d newValue2 = _mm256_fmadd_pd(tentative2,newAlpha2,oldDj);
3737        oldDj = _mm256_fmadd_pd(newAlpha2, tentative2, oldDj);
3738        __v4df bitsDj = _mm256_cmp_pd(oldDj, dualT2, _CMP_LT_OS);
3739        __v4df bitsAcceptable = _mm256_cmp_pd(newAlpha2, acceptable2, _CMP_GE_OS);
3740        _mm256_store_pd(reinterpret_cast< double * >(goodDj), bitsDj);
3741        _mm256_store_pd(reinterpret_cast< double * >(acceptableY), bitsAcceptable);
3742        _mm256_store_pd(alphaY, newAlpha2);
3743#ifndef USE_USE_AVX
3744#undef NDEBUG
3745        for (int i = 0; i < CHECK_CHUNK; i++) {
3746          assert(newValueX[i] > 0.0 == (goodDj[i]));
3747          //assert(acceptableX[i]==(acceptableY[i]));
3748          assert(oldValueX[i] == oldValueY[i]);
3749          assert(alphaX[i] == alphaY[i]);
3750        }
3751        for (int i = 0; i < CHECK_CHUNK; i++) {
3752          bool g1 = newValueX[i] < 0.0;
3753          bool g2 = goodDj[i] != 0;
3754          if (g1 != g2)
3755            abort();
3756          //if(acceptableX[i]!=(acceptableY[i]))abort();
3757          if (fabs(oldValueX[i] - oldValueY[i]) > 1.0e-5 + +(1.0e-10 * fabs(oldValueX[i])))
3758            abort();
3759          if (alphaX[i] != alphaY[i])
3760            abort();
3761        }
3762#endif
3763#endif
3764        for (int i = 0; i < CHECK_CHUNK + 1; i++) {
3765#ifndef USE_USE_AVX
3766          double newValue = newValueX[i];
3767          bool newStore = newValue < 0.0;
3768          if (store) {
3769            // add to list
3770            bool acceptable = acceptableX[i - 1];
3771            spare[numberRemaining] = work[i - 1];
3772            index[numberRemaining++] = which[i - 1] + addSequence;
3773            double value = oldValueX[i - 1] - upperTheta * alphaX[i - 1];
3774            if (value < dualT && acceptable) {
3775              upperTheta = (oldValueX[i - 1] - dualT) / alphaX[i - 1];
3776            }
3777          }
3778#else
3779          bool newStore = goodDj[i] != 0;
3780          if (store) {
3781            // add to list
3782            bool acceptable = acceptableY[i - 1];
3783            spare[numberRemaining] = work[i - 1];
3784            index[numberRemaining++] = which[i - 1] + addSequence;
3785            double value = oldValueY[i - 1] - upperTheta * alphaY[i - 1];
3786            if (value < dualT && acceptable) {
3787              upperTheta = (oldValueY[i - 1] - dualT) / alphaY[i - 1];
3788            }
3789          }
3790#endif
3791          store = newStore;
3792        }
3793        which += CHECK_CHUNK;
3794        work += CHECK_CHUNK;
3795      }
3796#endif
3797    }
3798#if ABOCA_LITE
3799    if (numberThreads) {
3800      work = columnArray->denseVector();
3801      number = columnArray->getNumElements();
3802      which = columnArray->getIndices();
3803      reducedCost = reducedCostWork_;
3804      unsigned char *statusArray = status_;
3805
3806      clpTempInfo info[ABOCA_LITE];
3807      int chunk = (number + numberThreads - 1) / numberThreads;
3808      int n = 0;
3809      int nR = numberRemaining;
3810      for (int i = 0; i < numberThreads; i++) {
3811        info[i].which = const_cast< int * >(which + n);
3812        info[i].work = const_cast< double * >(work + n);
3813        info[i].numberToDo = CoinMin(chunk, number - n);
3814        n += chunk;
3815        info[i].index = index + nR;
3816        info[i].spare = spare + nR;
3817        nR += chunk;
3818        info[i].reducedCost = const_cast< double * >(reducedCost);
3819        info[i].upperTheta = upperTheta;
3820        info[i].acceptablePivot = acceptablePivot;
3821        info[i].status = statusArray;
3822        info[i].tolerance = dualTolerance_;
3823      }
3824      // for gcc - get cilk out of function to stop avx2 error
3825      dualColumn000(numberThreads, info);
3826      moveAndZero(info, 1, NULL);
3827      for (int i = 0; i < numberThreads; i++) {
3828        numberRemaining += info[i].numberRemaining;
3829        upperTheta = CoinMin(upperTheta, static_cast< double >(info[i].upperTheta));
3830      }
3831    }
3832#endif
3833  } else {
3834    // some free or super basic
3835    for (int iSection = 0; iSection < 2; iSection++) {
3836
3837      int addSequence;
3838
3839      if (!iSection) {
3840        work = rowArray->denseVector();
3841        number = rowArray->getNumElements();
3842        which = rowArray->getIndices();
3843        reducedCost = rowReducedCost_;
3844        addSequence = numberColumns_;
3845      } else {
3846        work = columnArray->denseVector();
3847        number = columnArray->getNumElements();
3848        which = columnArray->getIndices();
3849        reducedCost = reducedCostWork_;
3850        addSequence = 0;
3851      }
3852
3853      for (i = 0; i < number; i++) {
3854        int iSequence = which[i];
3855        double alpha;
3856        double oldValue;
3857        double value;
3858        bool keep;
3859
3860        switch (getStatus(iSequence + addSequence)) {
3861
3862        case basic:
3863        case ClpSimplex::isFixed:
3864          break;
3865        case isFree:
3866        case superBasic:
3867          alpha = work[i];
3868          oldValue = reducedCost[iSequence];
3869          // If free has to be very large - should come in via dualRow
3870          //if (getStatus(iSequence+addSequence)==isFree&&fabs(alpha)<1.0e-3)
3871          //break;
3872          if (oldValue > dualTolerance_) {
3873            keep = true;
3874          } else if (oldValue < -dualTolerance_) {
3875            keep = true;
3876          } else {
3877            if (fabs(alpha) > CoinMax(10.0 * acceptablePivot, 1.0e-5)) {
3878              keep = true;
3879            } else {
3880              keep = false;
3881              badFree = CoinMax(badFree, fabs(alpha));
3882            }
3883          }
3884          if (keep) {
3885            // free - choose largest
3886            if (fabs(alpha) > freePivot) {
3887              freePivot = fabs(alpha);
3888              sequenceIn_ = iSequence + addSequence;
3889              theta_ = oldValue / alpha;
3890              alpha_ = alpha;
3891            }
3892            // give fake bounds if possible
3893            int jSequence = iSequence + addSequence;
3894            if (2.0 * fabs(solution_[jSequence]) < dualBound_) {
3895              FakeBound bound = getFakeBound(jSequence);
3896              assert(bound == ClpSimplexDual::noFake);
3897              setFakeBound(jSequence, ClpSimplexDual::bothFake);
3898              numberFake_++;
3899              value = oldValue - tentativeTheta * alpha;
3900              if (value > dualTolerance_) {
3901                // pretend coming in from upper bound
3902                upper_[jSequence] = solution_[jSequence];
3903                lower_[jSequence] = upper_[jSequence] - dualBound_;
3904                setColumnStatus(jSequence, ClpSimplex::atUpperBound);
3905              } else {
3906                // pretend coming in from lower bound
3907                lower_[jSequence] = solution_[jSequence];
3908                upper_[jSequence] = lower_[jSequence] + dualBound_;
3909                setColumnStatus(jSequence, ClpSimplex::atLowerBound);
3910              }
3911            }
3912          }
3913          break;
3914        case atUpperBound:
3915          alpha = work[i];
3916          oldValue = reducedCost[iSequence];
3917          value = oldValue - tentativeTheta * alpha;
3918          //assert (oldValue<=dualTolerance_*1.0001);
3919          if (value > dualTolerance_) {
3920            value = oldValue - upperTheta * alpha;
3921            if (value > dualTolerance_ && -alpha >= acceptablePivot) {
3922              upperTheta = (oldValue - dualTolerance_) / alpha;
3923              //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
3924            }
3925            // add to list
3926            spare[numberRemaining] = alpha;
3927            index[numberRemaining++] = iSequence + addSequence;
3928          }
3929          break;
3930        case atLowerBound:
3931          alpha = work[i];
3932          oldValue = reducedCost[iSequence];
3933          value = oldValue - tentativeTheta * alpha;
3934          //assert (oldValue>=-dualTolerance_*1.0001);
3935          if (value < -dualTolerance_) {
3936            value = oldValue - upperTheta * alpha;
3937            if (value < -dualTolerance_ && alpha >= acceptablePivot) {
3938              upperTheta = (oldValue + dualTolerance_) / alpha;
3939              //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
3940            }
3941            // add to list
3942            spare[numberRemaining] = alpha;
3943            index[numberRemaining++] = iSequence + addSequence;
3944          }
3945          break;
3946        }
3947      }
3948    }
3949  }
3950  upperReturn = upperTheta;
3951  return numberRemaining;
3952}
3953/*
3954   Row array has row part of pivot row (as duals so sign may be switched)
3955   Column array has column part.
3956   This chooses pivot column.
3957   Spare array will be needed when we start getting clever.
3958   We will check for basic so spare array will never overflow.
3959   If necessary will modify costs
3960*/
3961double
3962ClpSimplexDual::dualColumn(CoinIndexedVector *rowArray,
3963  CoinIndexedVector *columnArray,
3964  CoinIndexedVector *spareArray,
3965  CoinIndexedVector *spareArray2,
3966  double acceptablePivot,
3967  CoinBigIndex * /*dubiousWeights*/)
3968{
3969  int numberPossiblySwapped = 0;
3970  int numberRemaining = 0;
3971
3972  double totalThru = 0.0; // for when variables flip
3973  //double saveAcceptable=acceptablePivot;
3974  //acceptablePivot=1.0e-9;
3975
3976  double bestEverPivot = acceptablePivot;
3977  int lastSequence = -1;
3978  double lastPivot = 0.0;
3979  double upperTheta;
3980  double newTolerance = dualTolerance_;
3981  //newTolerance = dualTolerance_+1.0e-6*dblParam_[ClpDualTolerance];
3982  // will we need to increase tolerance
3983  //bool thisIncrease = false;
3984  // If we think we need to modify costs (not if something from broad sweep)
3985  bool modifyCosts = false;
3986  // Increase in objective due to swapping bounds (may be negative)
3987  double increaseInObjective = 0.0;
3988
3989  // use spareArrays to put ones looked at in
3990  // we are going to flip flop between
3991  int iFlip = 0;
3992  // Possible list of pivots
3993  int interesting[2];
3994  // where possible swapped ones are
3995  int swapped[2];
3996  // for zeroing out arrays after
3997  int marker[2][2];
3998  // pivot elements
3999  double *array[2], *spare, *spare2;
4000  // indices
4001  int *indices[2], *index, *index2;
4002  spareArray2->clear();
4003  array[0] = spareArray->denseVector();
4004  indices[0] = spareArray->getIndices();
4005  spare = array[0];
4006  index = indices[0];
4007  array[1] = spareArray2->denseVector();
4008  indices[1] = spareArray2->getIndices();
4009  int i;
4010
4011  // initialize lists
4012  for (i = 0; i < 2; i++) {
4013    interesting[i] = 0;
4014    swapped[i] = numberColumns_;
4015    marker[i][0] = 0;
4016    marker[i][1] = numberColumns_;
4017  }
4018  /*
4019       First we get a list of possible pivots.  We can also see if the
4020       problem looks infeasible or whether we want to pivot in free variable.
4021       This may make objective go backwards but can only happen a finite
4022       number of times and I do want free variables basic.
4023
4024       Then we flip back and forth.  At the start of each iteration
4025       interesting[iFlip] should have possible candidates and swapped[iFlip]
4026       will have pivots if we decide to take a previous pivot.
4027       At end of each iteration interesting[1-iFlip] should have
4028       candidates if we go through this theta and swapped[1-iFlip]
4029       pivots if we don't go through.
4030
4031       At first we increase theta and see what happens.  We start
4032       theta at a reasonable guess.  If in right area then we do bit by bit.
4033
4034      */
4035
4036  // do first pass to get possibles
4037  upperTheta = 1.0e31;
4038  double bestPossible = 1.0;
4039  double badFree = 0.0;
4040  alpha_ = 0.0;
4041  if (spareIntArray_[0] >= 0) {
4042    numberRemaining = dualColumn0(rowArray, columnArray, spareArray,
4043      acceptablePivot, upperTheta, badFree);
4044  } else {
4045    // already done
4046    numberRemaining = spareArray->getNumElements();
4047    spareArray->setNumElements(0);
4048    upperTheta = spareDoubleArray_[0];
4049    if (spareIntArray_[0] == -1) {
4050      theta_ = spareDoubleArray_[2];
4051      alpha_ = spareDoubleArray_[3];
4052      sequenceIn_ = spareIntArray_[1];
4053    } else {
4054#if 0
4055#undef NDEBUG
4056               int n = numberRemaining;
4057               double u = upperTheta;
4058               upperTheta = 1.0e31;
4059               CoinIndexedVector temp(4000);
4060               numberRemaining = dualColumn0(rowArray, columnArray, &temp,
4061                                             acceptablePivot, upperTheta, badFree);
4062               assert (n == numberRemaining);
4063               double * spare = spareArray->denseVector();
4064               int * index = spareArray->getIndices();
4065               double * spareX = temp.denseVector();
4066               int * indexX = temp.getIndices();
4067               CoinSort_2(spare,spare+n,index);
4068               CoinSort_2(spareX,spareX+n,indexX);
4069               for (int i=0;i<n;i++) {
4070                 assert (index[i]==indexX[i]);
4071                 assert (fabs(spare[i]-spareX[i])<1.0e-6);
4072               }
4073               assert (fabs(u - upperTheta) < 1.0e-7);
4074#endif
4075    }
4076  }
4077  // switch off
4078  spareIntArray_[0] = 0;
4079  // We can also see if infeasible or pivoting on free
4080  double tentativeTheta = 1.0e25;
4081  interesting[0] = numberRemaining;
4082  marker[0][0] = numberRemaining;
4083
4084  if (!numberRemaining && sequenceIn_ < 0)
4085    return 0.0; // Looks infeasible
4086
4087    // If sum of bad small pivots too much
4088#define MORE_CAREFUL
4089#ifdef MORE_CAREFUL
4090  bool badSumPivots = false;
4091#endif
4092  if (sequenceIn_ >= 0) {
4093    // free variable - always choose
4094  } else {
4095
4096    theta_ = 1.0e50;
4097    // now flip flop between spare arrays until reasonable theta
4098    tentativeTheta = CoinMax(10.0 * upperTheta, 1.0e-7);
4099
4100    // loops increasing tentative theta until can't go through
4101
4102    while (tentativeTheta < 1.0e22) {
4103      double thruThis = 0.0;
4104
4105      double bestPivot = acceptablePivot;
4106      int bestSequence = -1;
4107
4108      numberPossiblySwapped = numberColumns_;
4109      numberRemaining = 0;
4110
4111      upperTheta = 1.0e50;
4112
4113      spare = array[iFlip];
4114      index = indices[iFlip];
4115      spare2 = array[1 - iFlip];
4116      index2 = indices[1 - iFlip];
4117
4118      // try 3 different ways
4119      // 1 bias increase by ones with slightly wrong djs
4120      // 2 bias by all
4121      // 3 bias by all - tolerance
4122#define TRYBIAS 3
4123
4124      double increaseInThis = 0.0; //objective increase in this loop
4125
4126      for (i = 0; i < interesting[iFlip]; i++) {
4127        int iSequence = index[i];
4128        double alpha = spare[i];
4129        double oldValue = dj_[iSequence];
4130        double value = oldValue - tentativeTheta * alpha;
4131
4132        if (alpha < 0.0) {
4133          //at upper bound
4134          if (value > newTolerance) {
4135            double range = upper_[iSequence] - lower_[iSequence];
4136            thruThis -= range * alpha;
4137#if TRYBIAS == 1
4138            if (oldValue > 0.0)
4139              increaseInThis -= oldValue * range;
4140#elif TRYBIAS == 2
4141            increaseInThis -= oldValue * range;
4142#else
4143            increaseInThis -= (oldValue + dualTolerance_) * range;
4144#endif
4145            // goes on swapped list (also means candidates if too many)
4146            spare2[--numberPossiblySwapped] = alpha;
4147            index2[numberPossiblySwapped] = iSequence;
4148            if (fabs(alpha) > bestPivot) {
4149              bestPivot = fabs(alpha);
4150              bestSequence = numberPossiblySwapped;
4151            }
4152          } else {
4153            value = oldValue - upperTheta * alpha;
4154            if (value > newTolerance && -alpha >= acceptablePivot)
4155              upperTheta = (oldValue - newTolerance) / alpha;
4156            spare2[numberRemaining] = alpha;
4157            index2[numberRemaining++] = iSequence;
4158          }
4159        } else {
4160          // at lower bound
4161          if (value < -newTolerance) {
4162            double range = upper_[iSequence] - lower_[iSequence];
4163            thruThis += range * alpha;
4164            //?? is this correct - and should we look at good ones
4165#if TRYBIAS == 1
4166            if (oldValue < 0.0)
4167              increaseInThis += oldValue * range;
4168#elif TRYBIAS == 2
4169            increaseInThis += oldValue * range;
4170#else
4171            increaseInThis += (oldValue - dualTolerance_) * range;
4172#endif
4173            // goes on swapped list (also means candidates if too many)
4174            spare2[--numberPossiblySwapped] = alpha;
4175            index2[numberPossiblySwapped] = iSequence;
4176            if (fabs(alpha) > bestPivot) {
4177              bestPivot = fabs(alpha);
4178              bestSequence = numberPossiblySwapped;
4179            }
4180          } else {
4181            value = oldValue - upperTheta * alpha;
4182            if (value < -newTolerance && alpha >= acceptablePivot)
4183              upperTheta = (oldValue + newTolerance) / alpha;
4184            spare2[numberRemaining] = alpha;
4185            index2[numberRemaining++] = iSequence;
4186          }
4187        }
4188      }
4189      swapped[1 - iFlip] = numberPossiblySwapped;
4190      interesting[1 - iFlip] = numberRemaining;
4191      marker[1 - iFlip][0] = CoinMax(marker[1 - iFlip][0], numberRemaining);
4192      marker[1 - iFlip][1] = CoinMin(marker[1 - iFlip][1], numberPossiblySwapped);
4193
4194      double check = fabs(totalThru + thruThis);
4195      // add a bit
4196      check += 1.0e-8 + 1.0e-10 * check;
4197      if (check >= fabs(dualOut_) || increaseInObjective + increaseInThis < 0.0) {
4198        // We should be pivoting in this batch
4199        // so compress down to this lot
4200        numberRemaining = 0;
4201        for (i = numberColumns_ - 1; i >= swapped[1 - iFlip]; i--) {
4202          spare[numberRemaining] = spare2[i];
4203          index[numberRemaining++] = index2[i];
4204        }
4205        interesting[iFlip] = numberRemaining;
4206        int iTry;
4207#define MAXTRY 100
4208        // first get ratio with tolerance
4209        for (iTry = 0; iTry < MAXTRY; iTry++) {
4210
4211          upperTheta = 1.0e50;
4212          numberPossiblySwapped = numberColumns_;
4213          numberRemaining = 0;
4214
4215          increaseInThis = 0.0; //objective increase in this loop
4216
4217          thruThis = 0.0;
4218
4219          spare = array[iFlip];
4220          index = indices[iFlip];
4221          spare2 = array[1 - iFlip];
4222          index2 = indices[1 - iFlip];
4223          for (i = 0; i < interesting[iFlip]; i++) {
4224            int iSequence = index[i];
4225            double alpha = spare[i];
4226            double oldValue = dj_[iSequence];
4227            double value = oldValue - upperTheta * alpha;
4228
4229            if (alpha < 0.0) {
4230              //at upper bound
4231              if (value > newTolerance) {
4232                if (-alpha >= acceptablePivot) {
4233                  upperTheta = (oldValue - newTolerance) / alpha;
4234                }
4235              }
4236            } else {
4237              // at lower bound
4238              if (value < -newTolerance) {
4239                if (alpha >= acceptablePivot) {
4240                  upperTheta = (oldValue + newTolerance) / alpha;
4241                }
4242              }
4243            }
4244          }
4245          bestPivot = acceptablePivot;
4246          sequenceIn_ = -1;
4247#ifdef DUBIOUS_WEIGHTS
4248          double bestWeight = COIN_DBL_MAX;
4249#endif
4250          double largestPivot = acceptablePivot;
4251          // now choose largest and sum all ones which will go through
4252          //printf("XX it %d number %d\n",numberIterations_,interesting[iFlip]);
4253          // Sum of bad small pivots
4254#ifdef MORE_CAREFUL
4255          double sumBadPivots = 0.0;
4256          badSumPivots = false;
4257#endif
4258          // Make sure upperTheta will work (-O2 and above gives problems)
4259          upperTheta *= 1.0000000001;
4260          for (i = 0; i < interesting[iFlip]; i++) {
4261            int iSequence = index[i];
4262            double alpha = spare[i];
4263            double value = dj_[iSequence] - upperTheta * alpha;
4264            double badDj = 0.0;
4265
4266            bool addToSwapped = false;
4267
4268            if (alpha < 0.0) {
4269              //at upper bound
4270              if (value >= 0.0) {
4271                addToSwapped = true;
4272#if TRYBIAS == 1
4273                badDj = -CoinMax(dj_[iSequence], 0.0);
4274#elif TRYBIAS == 2
4275                badDj = -dj_[iSequence];
4276#else
4277                badDj = -dj_[iSequence] - dualTolerance_;
4278#endif
4279              }
4280            } else {
4281              // at lower bound
4282              if (value <= 0.0) {
4283                addToSwapped = true;
4284#if TRYBIAS == 1
4285                badDj = CoinMin(dj_[iSequence], 0.0);
4286#elif TRYBIAS == 2
4287                badDj = dj_[iSequence];
4288#else
4289                badDj = dj_[iSequence] - dualTolerance_;
4290#endif
4291              }
4292            }
4293            if (!addToSwapped) {
4294              // add to list of remaining
4295              spare2[numberRemaining] = alpha;
4296              index2[numberRemaining++] = iSequence;
4297            } else {
4298              // add to list of swapped
4299              spare2[--numberPossiblySwapped] = alpha;
4300              index2[numberPossiblySwapped] = iSequence;
4301              // select if largest pivot
4302              bool take = false;
4303              double absAlpha = fabs(alpha);
4304#ifdef DUBIOUS_WEIGHTS
4305              // User could do anything to break ties here
4306              double weight;
4307              if (dubiousWeights)
4308                weight = dubiousWeights[iSequence];
4309              else
4310                weight = 1.0;
4311              weight += randomNumberGenerator_.randomDouble() * 1.0e-2;
4312              if (absAlpha > 2.0 * bestPivot) {
4313                take = true;
4314              } else if (absAlpha > largestPivot) {
4315                // could multiply absAlpha and weight
4316                if (weight * bestPivot < bestWeight * absAlpha)
4317                  take = true;
4318              }
4319#else
4320              if (absAlpha > bestPivot)
4321                take = true;
4322#endif
4323#ifdef MORE_CAREFUL
4324              if (absAlpha < acceptablePivot && upperTheta < 1.0e20) {
4325                if (alpha < 0.0) {
4326                  //at upper bound
4327                  if (value > dualTolerance_) {
4328                    double gap = upper_[iSequence] - lower_[iSequence];
4329                    if (gap < 1.0e20)
4330                      sumBadPivots += value * gap;
4331                    else
4332                      sumBadPivots += 1.0e20;
4333                    //printf("bad %d alpha %g dj at upper %g\n",
4334                    //     iSequence,alpha,value);
4335                  }
4336                } else {
4337                  //at lower bound
4338                  if (value < -dualTolerance_) {
4339                    double gap = upper_[iSequence] - lower_[iSequence];
4340                    if (gap < 1.0e20)
4341                      sumBadPivots -= value * gap;
4342                    else
4343                      sumBadPivots += 1.0e20;
4344                    //printf("bad %d alpha %g dj at lower %g\n",
4345                    //     iSequence,alpha,value);
4346                  }
4347                }
4348              }
4349#endif
4350#ifdef FORCE_FOLLOW
4351              if (iSequence == force_in) {
4352                printf("taking %d - alpha %g best %g\n", force_in, absAlpha, largestPivot);
4353                take = true;
4354              }
4355#endif
4356              if (take) {
4357                sequenceIn_ = numberPossiblySwapped;
4358                bestPivot = absAlpha;
4359                theta_ = dj_[iSequence] / alpha;
4360                largestPivot = CoinMax(largestPivot, 0.5 * bestPivot);
4361#ifdef DUBIOUS_WEIGHTS
4362                bestWeight = weight;
4363#endif
4364                //printf(" taken seq %d alpha %g weight %d\n",
4365                //   iSequence,absAlpha,dubiousWeights[iSequence]);
4366              } else {
4367                //printf(" not taken seq %d alpha %g weight %d\n",
4368                //   iSequence,absAlpha,dubiousWeights[iSequence]);
4369              }
4370              double range = upper_[iSequence] - lower_[iSequence];
4371              thruThis += range * fabs(alpha);
4372              increaseInThis += badDj * range;
4373            }
4374          }
4375          marker[1 - iFlip][0] = CoinMax(marker[1 - iFlip][0], numberRemaining);
4376          marker[1 - iFlip][1] = CoinMin(marker[1 - iFlip][1], numberPossiblySwapped);
4377#ifdef MORE_CAREFUL
4378          // If we have done pivots and things look bad set alpha_ 0.0 to force factorization
4379          if (sumBadPivots > 1.0e4) {
4380            if (handler_->logLevel() > 1)
4381              *handler_ << "maybe forcing re-factorization - sum " << sumBadPivots << " " << factorization_->pivots() << " pivots" << CoinMessageEol;
4382            if (factorization_->pivots() > 3) {
4383              badSumPivots = true;
4384              break;
4385            }
4386          }
4387#endif
4388          swapped[1 - iFlip] = numberPossiblySwapped;
4389          interesting[1 - iFlip] = numberRemaining;
4390          // If we stop now this will be increase in objective (I think)
4391          double increase = (fabs(dualOut_) - totalThru) * theta_;
4392          increase += increaseInObjective;
4393          if (theta_ < 0.0)
4394            thruThis += fabs(dualOut_); // force using this one
4395          if (increaseInObjective < 0.0 && increase < 0.0 && lastSequence >= 0) {
4396            // back
4397            // We may need to be more careful - we could do by
4398            // switch so we always do fine grained?
4399            bestPivot = 0.0;
4400          } else {
4401            // add in
4402            totalThru += thruThis;
4403            increaseInObjective += increaseInThis;
4404          }
4405          if (bestPivot < 0.1 * bestEverPivot && bestEverPivot > 1.0e-6 && (bestPivot < 1.0e-3 || totalThru * 2.0 > fabs(dualOut_))) {
4406            // back to previous one
4407            sequenceIn_ = lastSequence;
4408            // swap regions
4409            iFlip = 1 - iFlip;
4410            break;
4411          } else if (sequenceIn_ == -1 && upperTheta > largeValue_) {
4412            if (lastPivot > acceptablePivot) {
4413              // back to previous one
4414              sequenceIn_ = lastSequence;
4415              // swap regions
4416              iFlip = 1 - iFlip;
4417            } else {
4418              // can only get here if all pivots too small
4419            }
4420            break;
4421          } else if (totalThru >= fabs(dualOut_)) {
4422            modifyCosts = true; // fine grain - we can modify costs
4423            break; // no point trying another loop
4424          } else {
4425            lastSequence = sequenceIn_;
4426            if (bestPivot > bestEverPivot)
4427              bestEverPivot = bestPivot;
4428            iFlip = 1 - iFlip;
4429            modifyCosts = true; // fine grain - we can modify costs
4430          }
4431        }
4432        if (iTry == MAXTRY)
4433          iFlip = 1 - iFlip; // flip back
4434        break;
4435      } else {
4436        // skip this lot
4437        if (bestPivot > 1.0e-3 || bestPivot > bestEverPivot) {
4438          bestEverPivot = bestPivot;
4439          lastSequence = bestSequence;
4440        } else {
4441          // keep old swapped
4442          CoinMemcpyN(array[iFlip] + swapped[iFlip],
4443            numberColumns_ - swapped[iFlip], array[1 - iFlip] + swapped[iFlip]);
4444          CoinMemcpyN(indices[iFlip] + swapped[iFlip],
4445            numberColumns_ - swapped[iFlip], indices[1 - iFlip] + swapped[iFlip]);
4446          marker[1 - iFlip][1] = CoinMin(marker[1 - iFlip][1], swapped[iFlip]);
4447          swapped[1 - iFlip] = swapped[iFlip];
4448        }
4449        increaseInObjective += increaseInThis;
4450        iFlip = 1 - iFlip; // swap regions
4451        tentativeTheta = 2.0 * upperTheta;
4452        totalThru += thruThis;
4453      }
4454    }
4455
4456    // can get here without sequenceIn_ set but with lastSequence
4457    if (sequenceIn_ < 0 && lastSequence >= 0) {
4458      // back to previous one
4459      sequenceIn_ = lastSequence;
4460      // swap regions
4461      iFlip = 1 - iFlip;
4462    }
4463
4464#define MINIMUMTHETA 1.0e-18
4465    // Movement should be minimum for anti-degeneracy - unless
4466    // fixed variable out
4467    double minimumTheta;
4468    if (upperOut_ > lowerOut_)
4469      minimumTheta = MINIMUMTHETA;
4470    else
4471      minimumTheta = 0.0;
4472    if (sequenceIn_ >= 0) {
4473      // at this stage sequenceIn_ is just pointer into index array
4474      // flip just so we can use iFlip
4475      iFlip = 1 - iFlip;
4476      spare = array[iFlip];
4477      index = indices[iFlip];
4478      double oldValue;
4479      alpha_ = spare[sequenceIn_];
4480      sequenceIn_ = indices[iFlip][sequenceIn_];
4481      oldValue = dj_[sequenceIn_];
4482      theta_ = CoinMax(oldValue / alpha_, 0.0);
4483      if (theta_ < minimumTheta && fabs(alpha_) < 1.0e5 && 1) {
4484        // can't pivot to zero
4485#if 0
4486                    if (oldValue - minimumTheta*alpha_ >= -dualTolerance_) {
4487                         theta_ = minimumTheta;
4488                    } else if (oldValue - minimumTheta*alpha_ >= -newTolerance) {
4489                         theta_ = minimumTheta;
4490                         thisIncrease = true;
4491                    } else {
4492                         theta_ = CoinMax((oldValue + newTolerance) / alpha_, 0.0);
4493                         thisIncrease = true;
4494                    }
4495#else
4496        theta_ = minimumTheta;
4497#endif
4498      }
4499      // may need to adjust costs so all dual feasible AND pivoted is exactly 0
4500      //int costOffset = numberRows_+numberColumns_;
4501      if (modifyCosts && !badSumPivots) {
4502        int i;
4503        for (i = numberColumns_ - 1; i >= swapped[iFlip]; i--) {
4504          int iSequence = index[i];
4505          double alpha = spare[i];
4506          double value = dj_[iSequence] - theta_ * alpha;
4507
4508          // can't be free here
4509
4510          if (alpha < 0.0) {
4511            //at upper bound
4512            if (value > dualTolerance_) {
4513              //thisIncrease = true;
4514#if CLP_CAN_HAVE_ZERO_OBJ < 2
4515#define MODIFYCOST 2
4516#endif
4517#if MODIFYCOST
4518              // modify cost to hit new tolerance
4519              double modification = alpha * theta_ - dj_[iSequence]
4520                + newTolerance;
4521              if ((specialOptions_ & (2048 + 4096 + 16384)) != 0) {
4522                if ((specialOptions_ & 16384) != 0) {
4523                  if (fabs(modification) < 1.0e-8)
4524                    modification = 0.0;
4525                } else if ((specialOptions_ & 2048) != 0) {
4526                  if (fabs(modification) < 1.0e-10)
4527                    modification = 0.0;
4528                } else {
4529                  if (fabs(modification) < 1.0e-12)
4530                    modification = 0.0;
4531                }
4532              }
4533              dj_[iSequence] += modification;
4534              cost_[iSequence] += modification;
4535              if (modification)
4536                numberChanged_++; // Say changed costs
4537                //cost_[iSequence+costOffset] += modification; // save change
4538#endif
4539            }
4540          } else {
4541            // at lower bound
4542            if (-value > dualTolerance_) {
4543              //thisIncrease = true;
4544#if MODIFYCOST
4545              // modify cost to hit new tolerance
4546              double modification = alpha * theta_ - dj_[iSequence]
4547                - newTolerance;
4548              //modification = CoinMax(modification,-dualTolerance_);
4549              //assert (fabs(modification)<1.0e-7);
4550              if ((specialOptions_ & (2048 + 4096)) != 0) {
4551                if ((specialOptions_ & 2048) != 0) {
4552                  if (fabs(modification) < 1.0e-10)
4553                    modification = 0.0;
4554                } else {
4555                  if (fabs(modification) < 1.0e-12)
4556                    modification = 0.0;
4557                }
4558              }
4559              dj_[iSequence] += modification;
4560              cost_[iSequence] += modification;
4561              if (modification)
4562                numberChanged_++; // Say changed costs
4563                //cost_[iSequence+costOffset] += modification; // save change
4564#endif
4565            }
4566          }
4567        }
4568      }
4569    }
4570  }
4571
4572#ifdef MORE_CAREFUL
4573  // If we have done pivots and things look bad set alpha_ 0.0 to force factorization
4574  if ((badSumPivots || fabs(theta_ * badFree) > 10.0 * dualTolerance_) && factorization_->pivots()) {
4575    if (handler_->logLevel() > 1)
4576      *handler_ << "forcing re-factorization" << CoinMessageEol;
4577    //printf("badSumPivots %g theta_ %g badFree %g\n",badSumPivots,theta_,badFree);
4578    sequenceIn_ = -1;
4579    acceptablePivot_ = -acceptablePivot_;
4580  }
4581#endif
4582  if (sequenceIn_ >= 0) {
4583    lowerIn_ = lower_[sequenceIn_];
4584    upperIn_ = upper_[sequenceIn_];
4585    valueIn_ = solution_[sequenceIn_];
4586    dualIn_ = dj_[sequenceIn_];
4587
4588    if (numberTimesOptimal_) {
4589      // can we adjust cost back closer to original
4590      //*** add coding
4591    }
4592#if MODIFYCOST > 1
4593    // modify cost to hit zero exactly
4594    // so (dualIn_+modification)==theta_*alpha_
4595    double modification = theta_ * alpha_ - dualIn_;
4596    // But should not move objective too much ??
4597#define DONT_MOVE_OBJECTIVE
4598#ifdef DONT_MOVE_OBJECTIVE
4599    double moveObjective = fabs(modification * solution_[sequenceIn_]);
4600    double smallMove = CoinMax(fabs(objectiveValue_), 1.0e-3);
4601    if (moveObjective > smallMove) {
4602      if (handler_->logLevel() > 1)
4603        printf("would move objective by %g - original mod %g sol value %g\n", moveObjective,
4604          modification, solution_[sequenceIn_]);
4605      modification *= smallMove / moveObjective;
4606    }
4607#endif
4608    if (badSumPivots)
4609      modification = 0.0;
4610    if ((specialOptions_ & (2048 + 4096)) != 0) {
4611      if ((specialOptions_ & 16384) != 0) {
4612        // in fast dual
4613        if (fabs(modification) < 1.0e-7)
4614          modification = 0.0;
4615      } else if ((specialOptions_ & 2048) != 0) {
4616        if (fabs(modification) < 1.0e-10)
4617          modification = 0.0;
4618      } else {
4619        if (fabs(modification) < 1.0e-12)
4620          modification = 0.0;
4621      }
4622    }
4623    dualIn_ += modification;
4624    dj_[sequenceIn_] = dualIn_;
4625    cost_[sequenceIn_] += modification;
4626    if (modification)
4627      numberChanged_++; // Say changed costs
4628      //int costOffset = numberRows_+numberColumns_;
4629      //cost_[sequenceIn_+costOffset] += modification; // save change
4630      //assert (fabs(modification)<1.0e-6);
4631#ifdef CLP_DEBUG
4632    if ((handler_->logLevel() & 32) && fabs(modification) > 1.0e-15)
4633      printf("exact %d new cost %g, change %g\n", sequenceIn_,
4634        cost_[sequenceIn_], modification);
4635#endif
4636#endif
4637
4638    if (alpha_ < 0.0) {
4639      // as if from upper bound
4640      directionIn_ = -1;
4641      upperIn_ = valueIn_;
4642    } else {
4643      // as if from lower bound
4644      directionIn_ = 1;
4645      lowerIn_ = valueIn_;
4646    }
4647    if (fabs(alpha_) < 1.0e-6) {
4648      // need bestPossible
4649      const double *work;
4650      int number;
4651      const int *which;
4652      const double *reducedCost;
4653      double tentativeTheta = 1.0e15;
4654      double upperTheta = 1.0e31;
4655      bestPossible = 0.0;
4656      double multiplier[] = { -1.0, 1.0 };
4657      double dualT = -dualTolerance_;
4658      int nSections = 2;
4659      int addSequence;
4660      for (int iSection = 0; iSection < nSections; iSection++) {
4661        if (!iSection) {
4662          work = rowArray->denseVector();
4663          number = rowArray->getNumElements();
4664          which = rowArray->getIndices();
4665          reducedCost = rowReducedCost_;
4666          addSequence = numberColumns_;
4667        } else {
4668          work = columnArray->denseVector();
4669          number = columnArray->getNumElements();
4670          which = columnArray->getIndices();
4671          reducedCost = reducedCostWork_;
4672          addSequence = 0;
4673        }
4674        for (i = 0; i < number; i++) {
4675          int iSequence = which[i];
4676          double alpha;
4677          double oldValue;
4678          double value;
4679          double mult = 1.0;
4680          switch (getStatus(iSequence + addSequence)) {
4681
4682          case basic:
4683          case ClpSimplex::isFixed:
4684            break;
4685          case isFree:
4686          case superBasic:
4687            alpha = work[i];
4688            bestPossible = CoinMax(bestPossible, fabs(alpha));
4689            break;
4690          case atUpperBound:
4691            mult = -1.0;
4692          case atLowerBound:
4693            alpha = work[i] * mult;
4694            if (alpha > 0.0) {
4695              oldValue = reducedCost[iSequence] * mult;
4696              value = oldValue - tentativeTheta * alpha;
4697              if (value < dualT) {
4698                bestPossible = CoinMax(bestPossible, alpha);
4699              }
4700            }
4701            break;
4702          }
4703        }
4704      }
4705    } else {
4706      bestPossible = fabs(alpha_);
4707    }
4708  } else {
4709    // no pivot
4710    bestPossible = 0.0;
4711    alpha_ = 0.0;
4712  }
4713  //if (thisIncrease)
4714  //dualTolerance_+= 1.0e-6*dblParam_[ClpDualTolerance];
4715
4716  // clear arrays
4717
4718  for (i = 0; i < 2; i++) {
4719    CoinZeroN(array[i], marker[i][0]);
4720    CoinZeroN(array[i] + marker[i][1], numberColumns_ - marker[i][1]);
4721  }
4722  return bestPossible;
4723}
4724#ifdef CLP_ALL_ONE_FILE
4725#undef MAXTRY
4726#endif
4727/* Checks if tentative optimal actually means unbounded
4728   Returns -3 if not, 2 if is unbounded */
4729int ClpSimplexDual::checkUnbounded(CoinIndexedVector *ray,
4730  CoinIndexedVector *spare,
4731  double changeCost)
4732{
4733  int status = 2; // say unbounded
4734  factorization_->updateColumn(spare, ray);
4735  // get reduced cost
4736  int i;
4737  int number = ray->getNumElements();
4738  int *index = ray->getIndices();
4739  double *array = ray->denseVector();
4740  for (i = 0; i < number; i++) {
4741    int iRow = index[i];
4742    int iPivot = pivotVariable_[iRow];
4743    changeCost -= cost(iPivot) * array[iRow];
4744  }
4745  double way;
4746  if (changeCost > 0.0) {
4747    //try going down
4748    way = 1.0;
4749  } else if (changeCost < 0.0) {
4750    //try going up
4751    way = -1.0;
4752  } else {
4753#ifdef CLP_DEBUG
4754    printf("can't decide on up or down\n");
4755#endif
4756    way = 0.0;
4757    status = -3;
4758  }
4759  double movement = 1.0e10 * way; // some largish number
4760  double zeroTolerance = 1.0e-14 * dualBound_;
4761  for (i = 0; i < number; i++) {
4762    int iRow = index[i];
4763    int iPivot = pivotVariable_[iRow];
4764    double arrayValue = array[iRow];
4765    if (fabs(arrayValue) < zeroTolerance)
4766      arrayValue = 0.0;
4767    double newValue = solution(iPivot) + movement * arrayValue;
4768    if (newValue > upper(iPivot) + primalTolerance_ || newValue < lower(iPivot) - primalTolerance_)
4769      status = -3; // not unbounded
4770  }
4771  if (status == 2) {
4772    // create ray
4773    delete[] ray_;
4774    ray_ = new double[numberColumns_];
4775    CoinZeroN(ray_, numberColumns_);
4776    for (i = 0; i < number; i++) {
4777      int iRow = index[i];
4778      int iPivot = pivotVariable_[iRow];
4779      double arrayValue = array[iRow];
4780      if (iPivot < numberColumns_ && fabs(arrayValue) >= zeroTolerance)
4781        ray_[iPivot] = way * array[iRow];
4782    }
4783  }
4784  ray->clear();
4785  return status;
4786}
4787//static int count_status=0;
4788//static double obj_status=0.0;
4789//static int check_status=123456789;//41754;
4790//static int count_alpha=0;
4791/* Checks if finished.  Updates status */
4792void ClpSimplexDual::statusOfProblemInDual(int &lastCleaned, int type,
4793  double *givenDuals, ClpDataSave &saveData,
4794  int ifValuesPass)
4795{
4796#ifdef CLP_INVESTIGATE_SERIAL
4797  if (z_thinks > 0 && z_thinks < 2)
4798    z_thinks += 2;
4799#endif
4800  bool arraysNotCreated = (type == 0);
4801  // If lots of iterations then adjust costs if large ones
4802  if (numberIterations_ > 4 * (numberRows_ + numberColumns_) && objectiveScale_ == 1.0) {
4803    double largest = 0.0;
4804    for (int i = 0; i < numberRows_; i++) {
4805      int iColumn = pivotVariable_[i];
4806      largest = CoinMax(largest, fabs(cost_[iColumn]));
4807    }
4808    if (largest > 1.0e6) {
4809      objectiveScale_ = 1.0e6 / largest;
4810      for (int i = 0; i < numberRows_ + numberColumns_; i++)
4811        cost_[i] *= objectiveScale_;
4812    }
4813  }
4814  int numberPivots = factorization_->pivots();
4815  double realDualInfeasibilities = 0.0;
4816  if (type == 2) {
4817    if (alphaAccuracy_ != -1.0)
4818      alphaAccuracy_ = -2.0;
4819    // trouble - restore solution
4820    CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4821    CoinMemcpyN(savedSolution_ + numberColumns_,
4822      numberRows_, rowActivityWork_);
4823    CoinMemcpyN(savedSolution_,
4824      numberColumns_, columnActivityWork_);
4825    // restore extra stuff
4826    int dummy;
4827    matrix_->generalExpanded(this, 6, dummy);
4828    forceFactorization_ = 1; // a bit drastic but ..
4829    changeMade_++; // say something changed
4830    // get correct bounds on all variables
4831    resetFakeBounds(0);
4832  }
4833  int tentativeStatus = problemStatus_;
4834  double changeCost;
4835  bool unflagVariables = true;
4836  bool weightsSaved = false;
4837  bool weightsSaved2 = numberIterations_ && !numberPrimalInfeasibilities_;
4838  int dontFactorizePivots = dontFactorizePivots_;
4839  if (type == 3) {
4840    type = 1;
4841    dontFactorizePivots = 1;
4842  }
4843  if (alphaAccuracy_ < 0.0 || !numberPivots || alphaAccuracy_ > 1.0e4 || numberPivots > 20) {
4844    if (problemStatus_ > -3 || numberPivots > dontFactorizePivots) {
4845      // factorize
4846      // later on we will need to recover from singularities
4847      // also we could skip if first time
4848      // save dual weights
4849      dualRowPivot_->saveWeights(this, 1);
4850      weightsSaved = true;
4851      if (type) {
4852        // is factorization okay?
4853        if (internalFactorize(1)) {
4854          // no - restore previous basis
4855          unflagVariables = false;
4856          assert(type == 1);
4857          changeMade_++; // say something changed
4858          // Keep any flagged variables
4859          int i;
4860          for (i = 0; i < numberRows_ + numberColumns_; i++) {
4861            if (flagged(i))
4862              saveStatus_[i] |= 64; //say flagged
4863          }
4864          CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4865          CoinMemcpyN(savedSolution_ + numberColumns_,
4866            numberRows_, rowActivityWork_);
4867          CoinMemcpyN(savedSolution_,
4868            numberColumns_, columnActivityWork_);
4869          // restore extra stuff
4870          int dummy;
4871          matrix_->generalExpanded(this, 6, dummy);
4872          // get correct bounds on all variables
4873          resetFakeBounds(1);
4874          // need to reject something
4875          char x = isColumn(sequenceOut_) ? 'C' : 'R';
4876          handler_->message(CLP_SIMPLEX_FLAG, messages_)
4877            << x << sequenceWithin(sequenceOut_)
4878            << CoinMessageEol;
4879#ifdef COIN_DEVELOP
4880          printf("flag d\n");
4881#endif
4882          setFlagged(sequenceOut_);
4883          progress_.clearBadTimes();
4884
4885          // Go to safe
4886          // not here - as can make more singular
4887          //factorization_->pivotTolerance(0.99);
4888          forceFactorization_ = 1; // a bit drastic but ..
4889          type = 2;
4890          //assert (internalFactorize(1)==0);
4891          if (internalFactorize(1)) {
4892            CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4893            CoinMemcpyN(savedSolution_ + numberColumns_,
4894              numberRows_, rowActivityWork_);
4895            CoinMemcpyN(savedSolution_,
4896              numberColumns_, columnActivityWork_);
4897            // restore extra stuff
4898            int dummy;
4899            matrix_->generalExpanded(this, 6, dummy);
4900            // debug
4901            int returnCode = internalFactorize(1);
4902            while (returnCode) {
4903              // ouch
4904              // switch off dense
4905              int saveDense = factorization_->denseThreshold();
4906              factorization_->setDenseThreshold(0);
4907              // Go to safe
4908              factorization_->pivotTolerance(0.99);
4909              // make sure will do safe factorization
4910              pivotVariable_[0] = -1;
4911              returnCode = internalFactorize(2);
4912              factorization_->setDenseThreshold(saveDense);
4913            }
4914            // get correct bounds on all variables
4915            resetFakeBounds(1);
4916          }
4917        }
4918      }
4919      if (problemStatus_ != -4 || numberPivots > 10)
4920        problemStatus_ = -3;
4921    }
4922  } else {
4923    //printf("testing with accuracy of %g and status of %d\n",alphaAccuracy_,problemStatus_);
4924    //count_alpha++;
4925    //if ((count_alpha%5000)==0)
4926    //printf("count alpha %d\n",count_alpha);
4927  }
4928  if (progress_.infeasibility_[0] < 1.0e-1 && primalTolerance_ == 1.0e-7 && progress_.iterationNumber_[0] > 0 && progress_.iterationNumber_[CLP_PROGRESS - 1] - progress_.iterationNumber_[0] > 25) {
4929    // default - so user did not set
4930    int iP;
4931    double minAverage = COIN_DBL_MAX;
4932    double maxAverage = 0.0;
4933    for (iP = 0; iP < CLP_PROGRESS; iP++) {
4934      int n = progress_.numberInfeasibilities_[iP];
4935      if (!n) {
4936        break;
4937      } else {
4938        double average = progress_.infeasibility_[iP];
4939        if (average > 0.1)
4940          break;
4941        average /= static_cast< double >(n);
4942        minAverage = CoinMin(minAverage, average);
4943        maxAverage = CoinMax(maxAverage, average);
4944      }
4945    }
4946    if (iP == CLP_PROGRESS && minAverage < 1.0e-5 && maxAverage < 1.0e-3) {
4947      // change tolerance
4948#if CBC_USEFUL_PRINTING > 0
4949      printf("CCchanging tolerance\n");
4950#endif
4951      primalTolerance_ = 1.0e-6;
4952      minimumPrimalTolerance_ = primalTolerance_;
4953      dblParam_[ClpPrimalTolerance] = 1.0e-6;
4954      moreSpecialOptions_ |= 4194304;
4955    }
4956  }
4957  // at this stage status is -3 or -4 if looks infeasible
4958  // get primal and dual solutions
4959#if 0
4960     {
4961          int numberTotal = numberRows_ + numberColumns_;
4962          double * saveSol = CoinCopyOfArray(solution_, numberTotal);
4963          double * saveDj = CoinCopyOfArray(dj_, numberTotal);
4964          double tolerance = type ? 1.0e-4 : 1.0e-8;
4965          // always if values pass
4966          double saveObj = objectiveValue_;
4967          double sumPrimal = sumPrimalInfeasibilities_;
4968          int numberPrimal = numberPrimalInfeasibilities_;
4969          double sumDual = sumDualInfeasibilities_;
4970          int numberDual = numberDualInfeasibilities_;
4971          gutsOfSolution(givenDuals, NULL);
4972          int j;
4973          double largestPrimal = tolerance;
4974          int iPrimal = -1;
4975          for (j = 0; j < numberTotal; j++) {
4976               double difference = solution_[j] - saveSol[j];
4977               if (fabs(difference) > largestPrimal) {
4978                    iPrimal = j;
4979                    largestPrimal = fabs(difference);
4980               }
4981          }
4982          double largestDual = tolerance;
4983          int iDual = -1;
4984          for (j = 0; j < numberTotal; j++) {
4985               double difference = dj_[j] - saveDj[j];
4986               if (fabs(difference) > largestDual && upper_[j] > lower_[j]) {
4987                    iDual = j;
4988                    largestDual = fabs(difference);
4989               }
4990          }
4991          if (!type) {
4992               if (fabs(saveObj - objectiveValue_) > 1.0e-5 ||
4993                         numberPrimal != numberPrimalInfeasibilities_ || numberPrimal != 1 ||
4994                         fabs(sumPrimal - sumPrimalInfeasibilities_) > 1.0e-5 || iPrimal >= 0 ||
4995                         numberDual != numberDualInfeasibilities_ || numberDual != 0 ||
4996                         fabs(sumDual - sumDualInfeasibilities_) > 1.0e-5 || iDual >= 0)
4997                    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",
4998                           type, numberIterations_, numberPivots,
4999                           numberPrimal, numberPrimalInfeasibilities_, sumPrimal, sumPrimalInfeasibilities_,
5000                           largestPrimal, iPrimal,
5001                           numberDual, numberDualInfeasibilities_, sumDual, sumDualInfeasibilities_,
5002                           largestDual, iDual,
5003                           saveObj, objectiveValue_);
5004          } else {
5005               if (fabs(saveObj - objectiveValue_) > 1.0e-5 ||
5006                         numberPrimalInfeasibilities_ || iPrimal >= 0 ||
5007                         numberDualInfeasibilities_ || iDual >= 0)
5008                    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",
5009                           type, numberIterations_, numberPivots,
5010                           numberPrimal, numberPrimalInfeasibilities_, sumPrimal, sumPrimalInfeasibilities_,
5011                           largestPrimal, iPrimal,
5012                           numberDual, numberDualInfeasibilities_, sumDual, sumDualInfeasibilities_,
5013                           largestDual, iDual,
5014                           saveObj, objectiveValue_);
5015          }
5016          delete [] saveSol;
5017          delete [] saveDj;
5018     }
5019#else
5020  if (type || ifValuesPass)
5021    gutsOfSolution(givenDuals, NULL);
5022#endif
5023  // If bad accuracy treat as singular
5024  if ((largestPrimalError_ > 1.0e15 || largestDualError_ > 1.0e15) && numberIterations_) {
5025    // restore previous basis
5026    unflagVariables = false;
5027    changeMade_++; // say something changed
5028    // Keep any flagged variables
5029    int i;
5030    for (i = 0; i < numberRows_ + numberColumns_; i++) {
5031      if (flagged(i))
5032        saveStatus_[i] |= 64; //say flagged
5033    }
5034    CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
5035    CoinMemcpyN(savedSolution_ + numberColumns_,
5036      numberRows_, rowActivityWork_);
5037    CoinMemcpyN(savedSolution_,
5038      numberColumns_, columnActivityWork_);
5039    // restore extra stuff
5040    int dummy;
5041    matrix_->generalExpanded(this, 6, dummy);
5042    // get correct bounds on all variables
5043    resetFakeBounds(1);
5044    // need to reject something
5045    int rejectedVariable = sequenceOut_;
5046    if (flagged(rejectedVariable)) {
5047      rejectedVariable = -1;
5048      for (int i = 0; i < numberRows_; i++) {
5049        int iSequence = pivotVariable_[i];
5050        if (!flagged(iSequence)) {
5051          rejectedVariable = iSequence;
5052          break;
5053        }
5054      }
5055      if (rejectedVariable < 0) {
5056        if (handler_->logLevel() > 1)
5057          printf("real trouble at line %d of ClpSimplexDual\n",
5058            __LINE__);
5059        problemStatus_ = 10;
5060        return;
5061      }
5062    }
5063    char x = isColumn(rejectedVariable) ? 'C' : 'R';
5064    handler_->message(CLP_SIMPLEX_FLAG, messages_)
5065      << x << sequenceWithin(rejectedVariable)
5066      << CoinMessageEol;
5067#ifdef COIN_DEVELOP
5068    printf("flag e\n");
5069#endif
5070    setFlagged(rejectedVariable);
5071    progress_.clearBadTimes();
5072
5073    // Go to safer
5074    double newTolerance = CoinMin(1.1 * factorization_->pivotTolerance(), 0.99);
5075    factorization_->pivotTolerance(newTolerance);
5076    forceFactorization_ = 1; // a bit drastic but ..
5077    if (alphaAccuracy_ != -1.0)
5078      alphaAccuracy_ = -2.0;
5079    type = 2;
5080    //assert (internalFactorize(1)==0);
5081    if (internalFactorize(1)) {
5082      CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
5083      CoinMemcpyN(savedSolution_ + numberColumns_,
5084        numberRows_, rowActivityWork_);
5085      CoinMemcpyN(savedSolution_,
5086        numberColumns_, columnActivityWork_);
5087      // restore extra stuff
5088      int dummy;
5089      matrix_->generalExpanded(this, 6, dummy);
5090      // debug
5091      int returnCode = internalFactorize(1);
5092      while (returnCode) {
5093        // ouch
5094        // switch off dense
5095        int saveDense = factorization_->denseThreshold();
5096        factorization_->setDenseThreshold(0);
5097        // Go to safe
5098        factorization_->pivotTolerance(0.99);
5099        // make sure will do safe factorization
5100        pivotVariable_[0] = -1;
5101        returnCode = internalFactorize(2);
5102        factorization_->setDenseThreshold(saveDense);
5103      }
5104      // get correct bounds on all variables
5105      resetFakeBounds(1);
5106    }
5107    // get primal and dual solutions
5108    gutsOfSolution(givenDuals, NULL);
5109  } else if (goodAccuracy()) {
5110    // Can reduce tolerance
5111    double newTolerance = CoinMax(0.995 * factorization_->pivotTolerance(), saveData.pivotTolerance_);
5112    factorization_->pivotTolerance(newTolerance);
5113  }
5114  bestObjectiveValue_ = CoinMax(bestObjectiveValue_,
5115    objectiveValue_ - bestPossibleImprovement_);
5116  bool reallyBadProblems = false;
5117  // Double check infeasibility if no action
5118  if (progress_.lastIterationNumber(0) == numberIterations_) {
5119    if (dualRowPivot_->looksOptimal()) {
5120      numberPrimalInfeasibilities_ = 0;
5121      sumPrimalInfeasibilities_ = 0.0;
5122    }
5123#if 1
5124  } else {
5125    double thisObj = objectiveValue_ - bestPossibleImprovement_;
5126#ifdef CLP_INVESTIGATE
5127    assert(bestPossibleImprovement_ > -1000.0 && objectiveValue_ > -1.0e100);
5128    if (bestPossibleImprovement_)
5129      printf("obj %g add in %g -> %g\n", objectiveValue_, bestPossibleImprovement_,
5130        thisObj);
5131#endif
5132    double lastObj = progress_.lastObjective(0);
5133#ifndef NDEBUG
5134#ifdef COIN_DEVELOP
5135    resetFakeBounds(-1);
5136#endif
5137#endif
5138#ifdef CLP_REPORT_PROGRESS
5139    ixxxxxx++;
5140    if (ixxxxxx >= ixxyyyy - 4 && ixxxxxx <= ixxyyyy) {
5141      char temp[20];
5142      sprintf(temp, "sol%d.out", ixxxxxx);
5143      printf("sol%d.out\n", ixxxxxx);
5144      FILE *fp = fopen(temp, "w");
5145      int nTotal = numberRows_ + numberColumns_;
5146      for (int i = 0; i < nTotal; i++)
5147        fprintf(fp, "%d %d %g %g %g %g %g\n",
5148          i, status_[i], lower_[i], solution_[i], upper_[i], cost_[i], dj_[i]);
5149      fclose(fp);
5150    }
5151#endif
5152    if (!ifValuesPass && firstFree_ < 0) {
5153      double testTol = 5.0e-3;
5154      if (progress_.timesFlagged() > 10) {
5155        testTol *= pow(2.0, progress_.timesFlagged() - 8);
5156      } else if (progress_.timesFlagged() > 5) {
5157        testTol *= 5.0;
5158      }
5159      if (lastObj > thisObj + testTol * (fabs(thisObj) + fabs(lastObj)) + testTol) {
5160        int maxFactor = factorization_->maximumPivots();
5161        if ((specialOptions_ & 1048576) == 0) {
5162          if (progress_.timesFlagged() > 10)
5163            progress_.incrementReallyBadTimes();
5164          if (maxFactor > 10 - 9) {
5165#ifdef COIN_DEVELOP
5166            printf("lastobj %g thisobj %g\n", lastObj, thisObj);
5167#endif
5168            //if (forceFactorization_<0)
5169            //forceFactorization_= maxFactor;
5170            //forceFactorization_ = CoinMax(1,(forceFactorization_>>1));
5171            if ((progressFlag_ & 4) == 0 && lastObj < thisObj + 1.0e4 && largestPrimalError_ < 1.0e2) {
5172              // Just save costs
5173              // save extra copy of cost_
5174              int nTotal = numberRows_ + numberColumns_;
5175              double *temp = new double[2 * nTotal];
5176              memcpy(temp, cost_, nTotal * sizeof(double));
5177              memcpy(temp + nTotal, cost_, nTotal * sizeof(double));
5178              delete[] cost_;
5179              cost_ = temp;
5180              objectiveWork_ = cost_;
5181              rowObjectiveWork_ = cost_ + numberColumns_;
5182              progressFlag_ |= 4;
5183            } else {
5184              forceFactorization_ = 1;
5185#ifdef COIN_DEVELOP
5186              printf("Reducing factorization frequency - bad backwards\n");
5187#endif
5188#if 1
5189              unflagVariables = false;
5190              changeMade_++; // say something changed
5191              int nTotal = numberRows_ + numberColumns_;
5192              CoinMemcpyN(saveStatus_, nTotal, status_);
5193              CoinMemcpyN(savedSolution_ + numberColumns_,
5194                numberRows_, rowActivityWork_);
5195              CoinMemcpyN(savedSolution_,
5196                numberColumns_, columnActivityWork_);
5197              if ((progressFlag_ & 4) == 0) {
5198                // save extra copy of cost_
5199                double *temp = new double[2 * nTotal];
5200                memcpy(temp, cost_, nTotal * sizeof(double));
5201                memcpy(temp + nTotal, cost_, nTotal * sizeof(double));
5202                delete[] cost_;
5203                cost_ = temp;
5204                objectiveWork_ = cost_;
5205                rowObjectiveWork_ = cost_ + numberColumns_;
5206                progressFlag_ |= 4;
5207              } else {
5208                memcpy(cost_, cost_ + nTotal, nTotal * sizeof(double));
5209              }
5210              // restore extra stuff
5211              int dummy;
5212              matrix_->generalExpanded(this, 6, dummy);
5213              double pivotTolerance = factorization_->pivotTolerance();
5214              if (pivotTolerance < 0.2)
5215                factorization_->pivotTolerance(0.2);
5216              else if (progress_.timesFlagged() > 2)
5217                factorization_->pivotTolerance(CoinMin(pivotTolerance * 1.1, 0.99));
5218              if (alphaAccuracy_ != -1.0)
5219                alphaAccuracy_ = -2.0;
5220              if (internalFactorize(1)) {
5221                CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
5222                CoinMemcpyN(savedSolution_ + numberColumns_,
5223                  numberRows_, rowActivityWork_);
5224                CoinMemcpyN(savedSolution_,
5225                  numberColumns_, columnActivityWork_);
5226                // restore extra stuff
5227                int dummy;
5228                matrix_->generalExpanded(this, 6, dummy);
5229                // debug
5230                int returnCode = internalFactorize(1);
5231                while (returnCode) {
5232                  // ouch
5233                  // switch off dense
5234                  int saveDense = factorization_->denseThreshold();
5235                  factorization_->setDenseThreshold(0);
5236                  // Go to safe
5237                  factorization_->pivotTolerance(0.99);
5238                  // make sure will do safe factorization
5239                  pivotVariable_[0] = -1;
5240                  returnCode = internalFactorize(2);
5241                  factorization_->setDenseThreshold(saveDense);
5242                }
5243              }
5244              resetFakeBounds(0);
5245              type = 2; // so will restore weights
5246              // get primal and dual solutions
5247              gutsOfSolution(givenDuals, NULL);
5248              if (numberPivots < 2) {
5249                // need to reject something
5250                char x = isColumn(sequenceOut_) ? 'C' : 'R';
5251                handler_->message(CLP_SIMPLEX_FLAG, messages_)
5252                  << x << sequenceWithin(sequenceOut_)
5253                  << CoinMessageEol;
5254#ifdef COIN_DEVELOP
5255                printf("flag d\n");
5256#endif
5257                setFlagged(sequenceOut_);
5258                progress_.clearBadTimes();
5259                progress_.incrementTimesFlagged();
5260              }
5261              if (numberPivots < 10)
5262                reallyBadProblems = true;
5263#ifdef COIN_DEVELOP
5264              printf("obj now %g\n", objectiveValue_);
5265#endif
5266              progress_.modifyObjective(objectiveValue_
5267                - bestPossibleImprovement_);
5268#endif
5269            }
5270          }
5271        } else {
5272          // in fast dual give up
5273#ifdef COIN_DEVELOP
5274          printf("In fast dual?\n");
5275#endif
5276          problemStatus_ = 3;
5277        }
5278      } else if (lastObj < thisObj - 1.0e-5 * CoinMax(fabs(thisObj), fabs(lastObj)) - 1.0e-3) {
5279        numberTimesOptimal_ = 0;
5280      }
5281    }
5282#endif
5283  }
5284  // Up tolerance if looks a bit odd
5285  if (numberIterations_ > CoinMax(1000, numberRows_ >> 4) && (specialOptions_ & 64) != 0) {
5286    if (sumPrimalInfeasibilities_ && sumPrimalInfeasibilities_ < 1.0e5) {
5287      int backIteration = progress_.lastIterationNumber(CLP_PROGRESS - 1);
5288      if (backIteration > 0 && numberIterations_ - backIteration < 9 * CLP_PROGRESS) {
5289        if (factorization_->pivotTolerance() < 0.9) {
5290          // up tolerance
5291          factorization_->pivotTolerance(CoinMin(factorization_->pivotTolerance() * 1.05 + 0.02, 0.91));
5292          //printf("tol now %g\n",factorization_->pivotTolerance());
5293          progress_.clearIterationNumbers();
5294        }
5295      }
5296    }
5297  }
5298  // Check if looping
5299  int loop;
5300  if (!givenDuals && type != 2)
5301    loop = progress_.looping();
5302  else
5303    loop = -1;
5304  if (progress_.reallyBadTimes() > 10) {
5305    problemStatus_ = 10; // instead - try other algorithm
5306#if COIN_DEVELOP > 2
5307    printf("returning at %d\n", __LINE__);
5308#endif
5309  }
5310  int situationChanged = 0;
5311  if (loop >= 0) {
5312    problemStatus_ = loop; //exit if in loop
5313    if (!problemStatus_) {
5314      // declaring victory
5315      numberPrimalInfeasibilities_ = 0;
5316      sumPrimalInfeasibilities_ = 0.0;
5317    } else {
5318      problemStatus_ = 10; // instead - try other algorithm
5319#if COIN_DEVELOP > 2
5320      printf("returning at %d\n", __LINE__);
5321#endif
5322    }
5323    return;
5324  } else if (loop < -1) {
5325    // something may have changed
5326    gutsOfSolution(NULL, NULL);
5327    situationChanged = 1;
5328  }
5329  // really for free variables in
5330  if ((progressFlag_ & 2) != 0) {
5331    situationChanged = 2;
5332  }
5333  progressFlag_ &= (~3); //reset progress flag
5334  if ((progressFlag_ & 4) != 0) {
5335    // save copy of cost_
5336    int nTotal = numberRows_ + numberColumns_;
5337    memcpy(cost_ + nTotal, cost_, nTotal * sizeof(double));
5338  }
5339  /*if (!numberIterations_&&sumDualInfeasibilities_)
5340       printf("OBJ %g sumPinf %g sumDinf %g\n",
5341        objectiveValue(),sumPrimalInfeasibilities_,
5342        sumDualInfeasibilities_);*/
5343  // mark as having gone optimal if looks like it
5344  if (!numberPrimalInfeasibilities_ && !numberDualInfeasibilities_)
5345    progressFlag_ |= 8;
5346  if (handler_->detail(CLP_SIMPLEX_STATUS, messages_) < 100) {
5347    handler_->message(CLP_SIMPLEX_STATUS, messages_)
5348      << numberIterations_ << objectiveValue();
5349    handler_->printing(sumPrimalInfeasibilities_ > 0.0)
5350      << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
5351    handler_->printing(sumDualInfeasibilities_ > 0.0)
5352      << sumDualInfeasibilities_ << numberDualInfeasibilities_;
5353    handler_->printing(numberDualInfeasibilitiesWithoutFree_
5354      < numberDualInfeasibilities_)
5355      << numberDualInfeasibilitiesWithoutFree_;
5356    handler_->message() << CoinMessageEol;
5357  }
5358#if 0
5359     count_status++;
5360     if (!numberIterations_)
5361       obj_status=-1.0e30;
5362     if (objectiveValue()<obj_status-0.01) {
5363       printf("Backward obj at %d from %g to %g\n",
5364              count_status,obj_status,objectiveValue());
5365     }
5366     obj_status=objectiveValue();
5367     if (count_status>=check_status-1) {
5368       printf("Trouble ahead - count_status %d\n",count_status);
5369     }
5370#endif
5371#if 0
5372     printf("IT %d %g %g(%d) %g(%d)\n",
5373            numberIterations_, objectiveValue(),
5374            sumPrimalInfeasibilities_, numberPrimalInfeasibilities_,
5375            sumDualInfeasibilities_, numberDualInfeasibilities_);
5376#endif
5377  double approximateObjective = objectiveValue_;
5378#ifdef CLP_REPORT_PROGRESS
5379  if (ixxxxxx >= ixxyyyy - 4 && ixxxxxx <= ixxyyyy) {
5380    char temp[20];
5381    sprintf(temp, "x_sol%d.out", ixxxxxx);
5382    FILE *fp = fopen(temp, "w");
5383    int nTotal = numberRows_ + numberColumns_;
5384    for (int i = 0; i < nTotal; i++)
5385      fprintf(fp, "%d %d %g %g %g %g %g\n",
5386        i, status_[i], lower_[i], solution_[i], upper_[i], cost_[i], dj_[i]);
5387    fclose(fp);
5388    if (ixxxxxx == ixxyyyy)
5389      exit(6);
5390  }
5391#endif
5392  realDualInfeasibilities = sumDualInfeasibilities_;
5393  double saveTolerance = dualTolerance_;
5394  // If we need to carry on cleaning variables
5395  if (!numberPrimalInfeasibilities_ && (specialOptions_ & 1024) != 0 && CLEAN_FIXED) {
5396    for (int iRow = 0; iRow < numberRows_; iRow++) {
5397      int iPivot = pivotVariable_[iRow];
5398      if (!flagged(iPivot) && pivoted(iPivot)) {
5399        // carry on
5400        numberPrimalInfeasibilities_ = -1;
5401        sumOfRelaxedPrimalInfeasibilities_ = 1.0;
5402        sumPrimalInfeasibilities_ = 1.0;
5403        break;
5404      }
5405    }
5406  }
5407  /* If we are primal feasible and any dual infeasibilities are on
5408        free variables then it is better to go to primal */
5409  if (!numberPrimalInfeasibilities_ && ((!numberDualInfeasibilitiesWithoutFree_ && numberDualInfeasibilities_) || (moreSpecialOptions_ & 16777216) != 0))
5410    problemStatus_ = 10;
5411  // dual bound coming in
5412  double saveDualBound = dualBound_;
5413  bool needCleanFake = false;
5414  while (problemStatus_ <= -3 && saveDualBound == dualBound_) {
5415    int cleanDuals = 0;
5416    if (situationChanged != 0)
5417      cleanDuals = 1;
5418    int numberChangedBounds = 0;
5419    int doOriginalTolerance = 0;
5420    if (lastCleaned == numberIterations_)
5421      doOriginalTolerance = 1;
5422    // check optimal
5423    // give code benefit of doubt
5424    if (sumOfRelaxedDualInfeasibilities_ == 0.0 && sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
5425      // say optimal (with these bounds etc)
5426      numberDualInfeasibilities_ = 0;
5427      sumDualInfeasibilities_ = 0.0;
5428      numberPrimalInfeasibilities_ = 0;
5429      sumPrimalInfeasibilities_ = 0.0;
5430    }
5431    //if (dualFeasible()||problemStatus_==-4||(primalFeasible()&&!numberDualInfeasibilitiesWithoutFree_)) {
5432    if (dualFeasible() || problemStatus_ == -4) {
5433      progress_.modifyObjective(objectiveValue_
5434        - bestPossibleImprovement_);
5435#ifdef COIN_DEVELOP
5436      if (sumDualInfeasibilities_ || bestPossibleImprovement_)
5437        printf("improve %g dualinf %g -> %g\n",
5438          bestPossibleImprovement_, sumDualInfeasibilities_,
5439          sumDualInfeasibilities_ * dualBound_);
5440#endif
5441      // see if cutoff reached
5442      double limit = 0.0;
5443      getDblParam(ClpDualObjectiveLimit, limit);
5444#if 0
5445               if(fabs(limit) < 1.0e30 && objectiveValue()*optimizationDirection_ >
5446                         limit + 1.0e-7 + 1.0e-8 * fabs(limit) && !numberAtFakeBound()) {
5447                    //looks infeasible on objective
5448                    if (perturbation_ == 101) {
5449                         cleanDuals = 1;
5450                         // Save costs
5451                         int numberTotal = numberRows_ + numberColumns_;
5452                         double * saveCost = CoinCopyOfArray(cost_, numberTotal);
5453                         // make sure fake bounds are back
5454                         changeBounds(1, NULL, changeCost);
5455                         createRim4(false);
5456                         // make sure duals are current
5457                         computeDuals(givenDuals);
5458                         checkDualSolution();
5459                         if(objectiveValue()*optimizationDirection_ >
5460                                   limit + 1.0e-7 + 1.0e-8 * fabs(limit) && !numberDualInfeasibilities_) {
5461                              perturbation_ = 102; // stop any perturbations
5462                              printf("cutoff test succeeded\n");
5463                         } else {
5464                              printf("cutoff test failed\n");
5465                              // put back
5466                              memcpy(cost_, saveCost, numberTotal * sizeof(double));
5467                              // make sure duals are current
5468                              computeDuals(givenDuals);
5469                              checkDualSolution();
5470                              progress_.modifyObjective(-COIN_DBL_MAX);
5471                              problemStatus_ = -1;
5472                         }
5473                         delete [] saveCost;
5474                    }
5475               }
5476#endif
5477      if (primalFeasible() && !givenDuals) {
5478        // may be optimal - or may be bounds are wrong
5479        handler_->message(CLP_DUAL_BOUNDS, messages_)
5480          << dualBound_
5481          << CoinMessageEol;
5482        // save solution in case unbounded
5483        double *saveColumnSolution = NULL;
5484        double *saveRowSolution = NULL;
5485        bool inCbc = (specialOptions_ & (0x01000000 | 16384)) != 0;
5486        if (!inCbc) {
5487          saveColumnSolution = CoinCopyOfArray(columnActivityWork_, numberColumns_);
5488          saveRowSolution = CoinCopyOfArray(rowActivityWork_, numberRows_);
5489        }
5490#ifndef COIN_MAX_DUAL_BOUND
5491#define COIN_MAX_DUAL_BOUND 1.0e20
5492#endif
5493        numberChangedBounds = (dualBound_ < COIN_MAX_DUAL_BOUND) ? changeBounds(0, rowArray_[3], changeCost) : 0;
5494        if (numberChangedBounds <= 0 && !numberDualInfeasibilities_) {
5495          //looks optimal - do we need to reset tolerance
5496          if (perturbation_ == 101) {
5497            perturbation_ = 102; // stop any perturbations
5498            cleanDuals = 1;
5499            // make sure fake bounds are back
5500            //computeObjectiveValue();
5501            changeBounds(1, NULL, changeCost);
5502            //computeObjectiveValue();
5503            createRim4(false);
5504            // make sure duals are current
5505            computeDuals(givenDuals);
5506            checkDualSolution();
5507            progress_.modifyObjective(-COIN_DBL_MAX);
5508#define DUAL_TRY_FASTER
5509#ifdef DUAL_TRY_FASTER
5510            if (numberDualInfeasibilities_) {
5511#endif
5512              numberChanged_ = 1; // force something to happen
5513              lastCleaned = numberIterations_ - 1;
5514#ifdef DUAL_TRY_FASTER
5515            } else {
5516              //double value = objectiveValue_;
5517              computeObjectiveValue(true);
5518              //printf("old %g new %g\n",value,objectiveValue_);
5519              //numberChanged_=1;
5520            }
5521#endif
5522          }
5523          if (lastCleaned < numberIterations_ && numberTimesOptimal_ < 4 && (numberChanged_ || (specialOptions_ & 4096) == 0)) {
5524#if CLP_CAN_HAVE_ZERO_OBJ
5525            if ((specialOptions_ & 16777216) == 0) {
5526#endif
5527              doOriginalTolerance = 2;
5528              numberTimesOptimal_++;
5529              changeMade_++; // say something changed
5530              if (numberTimesOptimal_ == 1) {
5531                dualTolerance_ = dblParam_[ClpDualTolerance];
5532              } else {
5533                if (numberTimesOptimal_ == 2) {
5534                  // better to have small tolerance even if slower
5535                  factorization_->zeroTolerance(CoinMin(factorization_->zeroTolerance(), 1.0e-15));
5536                }
5537                dualTolerance_ = dblParam_[ClpDualTolerance];
5538                dualTolerance_ *= pow(2.0, numberTimesOptimal_ - 1);
5539              }
5540              cleanDuals = 2; // If nothing changed optimal else primal
5541#if CLP_CAN_HAVE_ZERO_OBJ
5542            } else {
5543              // no cost - skip checks
5544              problemStatus_ = 0;
5545            }
5546#endif
5547          } else {
5548            problemStatus_ = 0; // optimal
5549            if (lastCleaned < numberIterations_ && numberChanged_) {
5550              handler_->message(CLP_SIMPLEX_GIVINGUP, messages_)
5551                << CoinMessageEol;
5552            }
5553          }
5554        } else {
5555          cleanDuals = 1;
5556          if (doOriginalTolerance == 1) {
5557            // check unbounded
5558            // find a variable with bad dj
5559            int iSequence;
5560            int iChosen = -1;
5561            if (!inCbc) {
5562              double largest = 100.0 * primalTolerance_;
5563              for (iSequence = 0; iSequence < numberRows_ + numberColumns_;
5564                   iSequence++) {
5565                double djValue = dj_[iSequence];
5566                double originalLo = originalLower(iSequence);
5567                double originalUp = originalUpper(iSequence);
5568                if (fabs(djValue) > fabs(largest)) {
5569                  if (getStatus(iSequence) != basic) {
5570                    if (djValue > 0 && originalLo < -1.0e20) {
5571                      if (djValue > fabs(largest)) {
5572                        largest = djValue;
5573                        iChosen = iSequence;
5574                      }
5575                    } else if (djValue < 0 && originalUp > 1.0e20) {
5576                      if (-djValue > fabs(largest)) {
5577                        largest = djValue;
5578                        iChosen = iSequence;
5579                      }
5580                    }
5581                  }
5582                }
5583              }
5584            }
5585            if (iChosen >= 0) {
5586              int iSave = sequenceIn_;
5587              sequenceIn_ = iChosen;
5588              unpack(rowArray_[1]);
5589              sequenceIn_ = iSave;
5590              // if dual infeasibilities then must be free vector so add in dual
5591              if (numberDualInfeasibilities_) {
5592                if (fabs(changeCost) > 1.0e-5)
5593                  COIN_DETAIL_PRINT(printf("Odd free/unbounded combo\n"));
5594                changeCost += cost_[iChosen];
5595              }
5596              problemStatus_ = checkUnbounded(rowArray_[1], rowArray_[0],
5597                changeCost);
5598              rowArray_[1]->clear();