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

Last change on this file since 2451 was 2451, checked in by forrest, 3 months ago

deal with memory leak and read error

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