# source:stable/1.17/Clp/src/ClpSimplexDual.cpp@2468

Last change on this file since 2468 was 2468, checked in by stefan, 15 months ago

sync with trunk r2467

• leave out MSVisualStudio/v15 for now: someone may still want to look into this; can wait until 1.18
• Property svn:eol-style set to `native`
• Property svn:keywords set to `Id`
File size: 295.9 KB
Line
1/* \$Id: ClpSimplexDual.cpp 2468 2019-05-03 04:31:14Z stefan \$ */
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
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
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
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
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)
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
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
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_);
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
1523#if ABOCA_LITE_FACTORIZATION
1525          cilk_sync;
1526        if (columnArray_[1]->getNumElements())
1528            pivotRow_, alpha_);
1530#endif
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)
1541        // if no pivots, bad update but reasonable alpha - take and invert
1542        if (updateStatus == 2 && !factorization_->pivots() && fabs(alpha_) > 1.0e-5)
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_);
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;
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) {
2113                  largest = diff - tol;
2114                }
2115              }
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
2204#endif
2206              int i;
2207              for (i = 0; i < numberRows_; i++) {
2208                if (rhs[i] < rowLowerWork_[i] - primalTolerance_ || rhs[i] > rowUpperWork_[i] + primalTolerance_) {
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
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_) {
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              }
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];
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
2436        clpTempInfo info[ABOCA_LITE];
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];
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];
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];
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];
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];
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];
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];
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];
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*/
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)) {
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);
2906        break;
2907      }
2908    } else {
2909      // use old
2910      forceThis = true;
2911    }
2912    if (!forceThis) {
2913      fclose(fpFollow);
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
2976#ifdef FORCE_FOLLOW
2977    if (forceThis) {
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_) {
3166              changeCost += movement * cost_[iSequence];
3167            } else {
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          }
3474          spare[numberRemaining] = alpha * mult;
3475          index[numberRemaining++] = iSequence;
3476        }
3477      }
3478    }
3479  }
3480  info.numberRemaining = numberRemaining;
3481  info.upperTheta = upperTheta;
3482}
3483static void
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{
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: {
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++) {
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++) {
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#elif defined(__arm__)
3560#include <arm_neon.h>
3561#else
3562#include <immintrin.h>
3563//#include <fmaintrin.h>
3564#endif
3565int ClpSimplexDual::dualColumn0(const CoinIndexedVector *rowArray,
3566  const CoinIndexedVector *columnArray,
3567  CoinIndexedVector *spareArray,
3568  double acceptablePivot,
3570{
3571  // do first pass to get possibles
3572  double *spare = spareArray->denseVector();
3573  int *index = spareArray->getIndices();
3574  const double *work;
3575  int number;
3576  const int *which;
3577  const double *reducedCost;
3578  // We can also see if infeasible or pivoting on free
3579  double tentativeTheta = 1.0e15;
3580  double upperTheta = 1.0e31;
3581  double freePivot = acceptablePivot;
3582  int numberRemaining = 0;
3583  int i;
3585  if ((moreSpecialOptions_ & 8) != 0) {
3586    // No free or super basic
3587    // bestPossible will re recomputed if necessary
3588#ifndef COIN_AVX2
3589    double multiplier[] = { -1.0, 1.0 };
3590#else
3591    double multiplier[4] = { 0.0, 0.0, -1.0, 1.0 };
3592#endif
3593    double dualT = -dualTolerance_;
3594#if ABOCA_LITE == 0
3595    int nSections = 2;
3596#else
3598    int nSections = numberThreads ? 1 : 2;
3599#endif
3600    for (int iSection = 0; iSection < nSections; iSection++) {
3601
3603      unsigned char *statusArray;
3604      if (!iSection) {
3605        work = rowArray->denseVector();
3606        number = rowArray->getNumElements();
3607        which = rowArray->getIndices();
3608        reducedCost = rowReducedCost_;
3610        statusArray = status_ + numberColumns_;
3611      } else {
3612        work = columnArray->denseVector();
3613        number = columnArray->getNumElements();
3614        which = columnArray->getIndices();
3615        reducedCost = reducedCostWork_;
3617        statusArray = status_;
3618      }
3619#ifndef COIN_AVX2
3620      for (i = 0; i < number; i++) {
3621        int iSequence = which[i];
3622        double alpha;
3623        double oldValue;
3624        double value;
3625
3626        assert(getStatus(iSequence + addSequence) != isFree
3627          && getStatus(iSequence + addSequence) != superBasic);
3628        int iStatus = (statusArray[iSequence] & 3) - 1;
3629        if (iStatus) {
3630          double mult = multiplier[iStatus - 1];
3631          alpha = work[i] * mult;
3632          if (alpha > 0.0) {
3633            oldValue = reducedCost[iSequence] * mult;
3634            value = oldValue - tentativeTheta * alpha;
3635            if (value < dualT) {
3636              value = oldValue - upperTheta * alpha;
3637              if (value < dualT && alpha >= acceptablePivot) {
3638                upperTheta = (oldValue - dualT) / alpha;
3639                //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
3640              }
3642              spare[numberRemaining] = alpha * mult;
3643              index[numberRemaining++] = iSequence + addSequence;
3644            }
3645          }
3646        }
3647      }
3648      //
3649#else
3650      //#define COIN_AVX2 4 // temp
3651#if COIN_AVX2 == 1
3652#define COIN_AVX2_SHIFT 0
3653#elif COIN_AVX2 == 2
3654#define COIN_AVX2_SHIFT 1
3655#elif COIN_AVX2 == 4
3656#define COIN_AVX2_SHIFT 2
3657#elif COIN_AVX2 == 8
3658#define COIN_AVX2_SHIFT 3
3659#else
3660      error;
3661#endif
3662      //#define COIN_ALIGN 8*COIN_AVX2 // later
3663      //#define COIN_ALIGN_DOUBLE COIN_AVX2
3664#define CHECK_CHUNK 4
3665      // round up
3666      int *whichX = const_cast< int * >(which);
3667      double *workX = const_cast< double * >(work);
3668      int nBlocks = (number + CHECK_CHUNK - 1) / CHECK_CHUNK;
3669      int n = nBlocks * CHECK_CHUNK + 1;
3670      for (int i = number; i < n; i++) {
3671        workX[i] = 0.0;
3672        whichX[i] = 0; // alpha will be zero so not chosen
3673      }
3674      bool acceptableX[CHECK_CHUNK + 1];
3675      double oldValueX[CHECK_CHUNK + 1];
3676      double newValueX[CHECK_CHUNK + 1];
3677      double alphaX[CHECK_CHUNK + 1];
3678      newValueX[CHECK_CHUNK] = 0.0;
3679#define USE_USE_AVX
3680      //#define CHECK_H 1
3681#ifdef USE_USE_AVX
3682#define NEED_AVX
3683#elif CHECK_H
3684#define NEED_AVX
3685#endif
3686#ifdef NEED_AVX
3687      double mult2[CHECK_CHUNK] __attribute__((aligned(64)));
3688      CoinInt64 acceptableY[CHECK_CHUNK] __attribute__((aligned(64)));
3689      CoinInt64 goodDj[CHECK_CHUNK + 1] __attribute__((aligned(64)));
3690      double oldValueY[CHECK_CHUNK] __attribute__((aligned(64)));
3691      double alphaY[CHECK_CHUNK + 1] __attribute__((aligned(64)));
3692      memset(acceptableY, 0, sizeof(acceptableY));
3693      memset(goodDj, 0, sizeof(goodDj));
3694      memset(oldValueY, 0, sizeof(oldValueY));
3695      memset(alphaY, 0, sizeof(alphaY));
3696      __m256d tentative2 = _mm256_set1_pd(-tentativeTheta);
3697      __m256d dualT2 = _mm256_set1_pd(dualT);
3698      __m256d acceptable2 = _mm256_set1_pd(acceptablePivot);
3699#endif
3700      for (int iBlock = 0; iBlock < nBlocks; iBlock++) {
3701        bool store = false;
3702        double alpha = 0.0;
3703        double oldValue = 0.0;
3704        double newValue = 0.0;
3705        double trueAlpha = 0.0;
3706        int jSequence = 0;
3707#ifndef USE_USE_AVX
3708        for (int i = 0; i < CHECK_CHUNK + 1; i++) {
3709          int iSequence = which[i];
3710          int iStatus = (statusArray[iSequence] & 3);
3711          double mult = multiplier[iStatus];
3712          double newAlpha = work[i] * mult;
3713          double oldDj = reducedCost[iSequence] * mult;
3714          newValue = (oldDj - tentativeTheta * newAlpha) - dualT;
3715          acceptableX[i] = newAlpha >= acceptablePivot;
3716          oldValueX[i] = oldDj;
3717          newValueX[i] = newValue;
3718          alphaX[i] = newAlpha;
3719        }
3720#endif
3721#ifdef NEED_AVX
3722        __m128i columns = _mm_load_si128((const __m128i *)which);
3723        // what do we get - this must be wrong
3724        // probably only 1 and 2 - can we be clever
3725        // fix
3726        //__m128i status; // = _mm256_i32gather_ps(statusArray,columns,1);
3727        //status.m128i_i32[0]=statusArray[columns.m128i_i32[0]];
3728        for (int i = 0; i < CHECK_CHUNK; i++) {
3729          int iSequence = which[i];
3730          int iStatus = (statusArray[iSequence] & 3);
3731          mult2[i] = multiplier[iStatus];
3732        }
3733        //__m256d newAlpha2 = _mm256_i32gather_pd(multiplier,status,1); // mult here
3735        __m256d oldDj = _mm256_i32gather_pd(reducedCost, columns, 8);
3736        oldDj = _mm256_mul_pd(oldDj, newAlpha2); // remember newAlpha==mult
3737        _mm256_store_pd(oldValueY, oldDj); // redo later
3739        newAlpha2 = _mm256_mul_pd(newAlpha2, work2); // now really newAlpha
3741        oldDj = _mm256_fmadd_pd(newAlpha2, tentative2, oldDj);
3742        __v4df bitsDj = _mm256_cmp_pd(oldDj, dualT2, _CMP_LT_OS);
3743        __v4df bitsAcceptable = _mm256_cmp_pd(newAlpha2, acceptable2, _CMP_GE_OS);
3744        _mm256_store_pd(reinterpret_cast< double * >(goodDj), bitsDj);
3745        _mm256_store_pd(reinterpret_cast< double * >(acceptableY), bitsAcceptable);
3746        _mm256_store_pd(alphaY, newAlpha2);
3747#ifndef USE_USE_AVX
3748#undef NDEBUG
3749        for (int i = 0; i < CHECK_CHUNK; i++) {
3750          assert(newValueX[i] > 0.0 == (goodDj[i]));
3751          //assert(acceptableX[i]==(acceptableY[i]));
3752          assert(oldValueX[i] == oldValueY[i]);
3753          assert(alphaX[i] == alphaY[i]);
3754        }
3755        for (int i = 0; i < CHECK_CHUNK; i++) {
3756          bool g1 = newValueX[i] < 0.0;
3757          bool g2 = goodDj[i] != 0;
3758          if (g1 != g2)
3759            abort();
3760          //if(acceptableX[i]!=(acceptableY[i]))abort();
3761          if (fabs(oldValueX[i] - oldValueY[i]) > 1.0e-5 + +(1.0e-10 * fabs(oldValueX[i])))
3762            abort();
3763          if (alphaX[i] != alphaY[i])
3764            abort();
3765        }
3766#endif
3767#endif
3768        for (int i = 0; i < CHECK_CHUNK + 1; i++) {
3769#ifndef USE_USE_AVX
3770          double newValue = newValueX[i];
3771          bool newStore = newValue < 0.0;
3772          if (store) {
3774            bool acceptable = acceptableX[i - 1];
3775            spare[numberRemaining] = work[i - 1];
3776            index[numberRemaining++] = which[i - 1] + addSequence;
3777            double value = oldValueX[i - 1] - upperTheta * alphaX[i - 1];
3778            if (value < dualT && acceptable) {
3779              upperTheta = (oldValueX[i - 1] - dualT) / alphaX[i - 1];
3780            }
3781          }
3782#else
3783          bool newStore = goodDj[i] != 0;
3784          if (store) {
3786            bool acceptable = acceptableY[i - 1];
3787            spare[numberRemaining] = work[i - 1];
3788            index[numberRemaining++] = which[i - 1] + addSequence;
3789            double value = oldValueY[i - 1] - upperTheta * alphaY[i - 1];
3790            if (value < dualT && acceptable) {
3791              upperTheta = (oldValueY[i - 1] - dualT) / alphaY[i - 1];
3792            }
3793          }
3794#endif
3795          store = newStore;
3796        }
3797        which += CHECK_CHUNK;
3798        work += CHECK_CHUNK;
3799      }
3800#endif
3801    }
3802#if ABOCA_LITE
3804      work = columnArray->denseVector();
3805      number = columnArray->getNumElements();
3806      which = columnArray->getIndices();
3807      reducedCost = reducedCostWork_;
3808      unsigned char *statusArray = status_;
3809
3810      clpTempInfo info[ABOCA_LITE];
3812      int n = 0;
3813      int nR = numberRemaining;
3814      for (int i = 0; i < numberThreads; i++) {
3815        info[i].which = const_cast< int * >(which + n);
3816        info[i].work = const_cast< double * >(work + n);
3817        info[i].numberToDo = CoinMin(chunk, number - n);
3818        n += chunk;
3819        info[i].index = index + nR;
3820        info[i].spare = spare + nR;
3821        nR += chunk;
3822        info[i].reducedCost = const_cast< double * >(reducedCost);
3823        info[i].upperTheta = upperTheta;
3824        info[i].acceptablePivot = acceptablePivot;
3825        info[i].status = statusArray;
3826        info[i].tolerance = dualTolerance_;
3827      }
3828      // for gcc - get cilk out of function to stop avx2 error
3830      moveAndZero(info, 1, NULL);
3831      for (int i = 0; i < numberThreads; i++) {
3832        numberRemaining += info[i].numberRemaining;
3833        upperTheta = CoinMin(upperTheta, static_cast< double >(info[i].upperTheta));
3834      }
3835    }
3836#endif
3837  } else {
3838    // some free or super basic
3839    for (int iSection = 0; iSection < 2; iSection++) {
3840
3842
3843      if (!iSection) {
3844        work = rowArray->denseVector();
3845        number = rowArray->getNumElements();
3846        which = rowArray->getIndices();
3847        reducedCost = rowReducedCost_;
3849      } else {
3850        work = columnArray->denseVector();
3851        number = columnArray->getNumElements();
3852        which = columnArray->getIndices();
3853        reducedCost = reducedCostWork_;
3855      }
3856
3857      for (i = 0; i < number; i++) {
3858        int iSequence = which[i];
3859        double alpha;
3860        double oldValue;
3861        double value;
3862        bool keep;
3863
3864        switch (getStatus(iSequence + addSequence)) {
3865
3866        case basic:
3867        case ClpSimplex::isFixed:
3868          break;
3869        case isFree:
3870        case superBasic:
3871          alpha = work[i];
3872          oldValue = reducedCost[iSequence];
3873          // If free has to be very large - should come in via dualRow
3875          //break;
3876          if (oldValue > dualTolerance_) {
3877            keep = true;
3878          } else if (oldValue < -dualTolerance_) {
3879            keep = true;
3880          } else {
3881            if (fabs(alpha) > CoinMax(10.0 * acceptablePivot, 1.0e-5)) {
3882              keep = true;
3883            } else {
3884              keep = false;
3886            }
3887          }
3888          if (keep) {
3889            // free - choose largest
3890            if (fabs(alpha) > freePivot) {
3891              freePivot = fabs(alpha);
3892              sequenceIn_ = iSequence + addSequence;
3893              theta_ = oldValue / alpha;
3894              alpha_ = alpha;
3895            }
3896            // give fake bounds if possible
3897            int jSequence = iSequence + addSequence;
3898            if (2.0 * fabs(solution_[jSequence]) < dualBound_) {
3899              FakeBound bound = getFakeBound(jSequence);
3900              assert(bound == ClpSimplexDual::noFake);
3901              setFakeBound(jSequence, ClpSimplexDual::bothFake);
3902              numberFake_++;
3903              value = oldValue - tentativeTheta * alpha;
3904              if (value > dualTolerance_) {
3905                // pretend coming in from upper bound
3906                upper_[jSequence] = solution_[jSequence];
3907                lower_[jSequence] = upper_[jSequence] - dualBound_;
3908                setColumnStatus(jSequence, ClpSimplex::atUpperBound);
3909              } else {
3910                // pretend coming in from lower bound
3911                lower_[jSequence] = solution_[jSequence];
3912                upper_[jSequence] = lower_[jSequence] + dualBound_;
3913                setColumnStatus(jSequence, ClpSimplex::atLowerBound);
3914              }
3915            }
3916          }
3917          break;
3918        case atUpperBound:
3919          alpha = work[i];
3920          oldValue = reducedCost[iSequence];
3921          value = oldValue - tentativeTheta * alpha;
3922          //assert (oldValue<=dualTolerance_*1.0001);
3923          if (value > dualTolerance_) {
3924            value = oldValue - upperTheta * alpha;
3925            if (value > dualTolerance_ && -alpha >= acceptablePivot) {
3926              upperTheta = (oldValue - dualTolerance_) / alpha;
3927              //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
3928            }
3930            spare[numberRemaining] = alpha;
3931            index[numberRemaining++] = iSequence + addSequence;
3932          }
3933          break;
3934        case atLowerBound:
3935          alpha = work[i];
3936          oldValue = reducedCost[iSequence];
3937          value = oldValue - tentativeTheta * alpha;
3938          //assert (oldValue>=-dualTolerance_*1.0001);
3939          if (value < -dualTolerance_) {
3940            value = oldValue - upperTheta * alpha;
3941            if (value < -dualTolerance_ && alpha >= acceptablePivot) {
3942              upperTheta = (oldValue + dualTolerance_) / alpha;
3943              //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
3944            }
3946            spare[numberRemaining] = alpha;
3947            index[numberRemaining++] = iSequence + addSequence;
3948          }
3949          break;
3950        }
3951      }
3952    }
3953  }
3954  upperReturn = upperTheta;
3955  return numberRemaining;
3956}
3957/*
3958   Row array has row part of pivot row (as duals so sign may be switched)
3959   Column array has column part.
3960   This chooses pivot column.
3961   Spare array will be needed when we start getting clever.
3962   We will check for basic so spare array will never overflow.
3963   If necessary will modify costs
3964*/
3965double
3966ClpSimplexDual::dualColumn(CoinIndexedVector *rowArray,
3967  CoinIndexedVector *columnArray,
3968  CoinIndexedVector *spareArray,
3969  CoinIndexedVector *spareArray2,
3970  double acceptablePivot,
3971  CoinBigIndex * /*dubiousWeights*/)
3972{
3973  int numberPossiblySwapped = 0;
3974  int numberRemaining = 0;
3975
3976  double totalThru = 0.0; // for when variables flip
3977  //double saveAcceptable=acceptablePivot;
3978  //acceptablePivot=1.0e-9;
3979
3980  double bestEverPivot = acceptablePivot;
3981  int lastSequence = -1;
3982  double lastPivot = 0.0;
3983  double upperTheta;
3984  double newTolerance = dualTolerance_;
3985  //newTolerance = dualTolerance_+1.0e-6*dblParam_[ClpDualTolerance];
3986  // will we need to increase tolerance
3987  //bool thisIncrease = false;
3988  // If we think we need to modify costs (not if something from broad sweep)
3989  bool modifyCosts = false;
3990  // Increase in objective due to swapping bounds (may be negative)
3991  double increaseInObjective = 0.0;
3992
3993  // use spareArrays to put ones looked at in
3994  // we are going to flip flop between
3995  int iFlip = 0;
3996  // Possible list of pivots
3997  int interesting[2];
3998  // where possible swapped ones are
3999  int swapped[2];
4000  // for zeroing out arrays after
4001  int marker[2][2];
4002  // pivot elements
4003  double *array[2], *spare, *spare2;
4004  // indices
4005  int *indices[2], *index, *index2;
4006  spareArray2->clear();
4007  array[0] = spareArray->denseVector();
4008  indices[0] = spareArray->getIndices();
4009  spare = array[0];
4010  index = indices[0];
4011  array[1] = spareArray2->denseVector();
4012  indices[1] = spareArray2->getIndices();
4013  int i;
4014
4015  // initialize lists
4016  for (i = 0; i < 2; i++) {
4017    interesting[i] = 0;
4018    swapped[i] = numberColumns_;
4019    marker[i][0] = 0;
4020    marker[i][1] = numberColumns_;
4021  }
4022  /*
4023       First we get a list of possible pivots.  We can also see if the
4024       problem looks infeasible or whether we want to pivot in free variable.
4025       This may make objective go backwards but can only happen a finite
4026       number of times and I do want free variables basic.
4027
4028       Then we flip back and forth.  At the start of each iteration
4029       interesting[iFlip] should have possible candidates and swapped[iFlip]
4030       will have pivots if we decide to take a previous pivot.
4031       At end of each iteration interesting[1-iFlip] should have
4032       candidates if we go through this theta and swapped[1-iFlip]
4033       pivots if we don't go through.
4034
4035       At first we increase theta and see what happens.  We start
4036       theta at a reasonable guess.  If in right area then we do bit by bit.
4037
4038      */
4039
4040  // do first pass to get possibles
4041  upperTheta = 1.0e31;
4042  double bestPossible = 1.0;
4044  alpha_ = 0.0;
4045  if (spareIntArray_[0] >= 0) {
4046    numberRemaining = dualColumn0(rowArray, columnArray, spareArray,
4048  } else {
4050    numberRemaining = spareArray->getNumElements();
4051    spareArray->setNumElements(0);
4052    upperTheta = spareDoubleArray_[0];
4053    if (spareIntArray_[0] == -1) {
4054      theta_ = spareDoubleArray_[2];
4055      alpha_ = spareDoubleArray_[3];
4056      sequenceIn_ = spareIntArray_[1];
4057    } else {
4058#if 0
4059#undef NDEBUG
4060               int n = numberRemaining;
4061               double u = upperTheta;
4062               upperTheta = 1.0e31;
4063               CoinIndexedVector temp(4000);
4064               numberRemaining = dualColumn0(rowArray, columnArray, &temp,
4066               assert (n == numberRemaining);
4067               double * spare = spareArray->denseVector();
4068               int * index = spareArray->getIndices();
4069               double * spareX = temp.denseVector();
4070               int * indexX = temp.getIndices();
4071               CoinSort_2(spare,spare+n,index);
4072               CoinSort_2(spareX,spareX+n,indexX);
4073               for (int i=0;i<n;i++) {
4074                 assert (index[i]==indexX[i]);
4075                 assert (fabs(spare[i]-spareX[i])<1.0e-6);
4076               }
4077               assert (fabs(u - upperTheta) < 1.0e-7);
4078#endif
4079    }
4080  }
4081  // switch off
4082  spareIntArray_[0] = 0;
4083  // We can also see if infeasible or pivoting on free
4084  double tentativeTheta = 1.0e25;
4085  interesting[0] = numberRemaining;
4086  marker[0][0] = numberRemaining;
4087
4088  if (!numberRemaining && sequenceIn_ < 0)
4089    return 0.0; // Looks infeasible
4090
4091    // If sum of bad small pivots too much
4092#define MORE_CAREFUL
4093#ifdef MORE_CAREFUL
4095#endif
4096  if (sequenceIn_ >= 0) {
4097    // free variable - always choose
4098  } else {
4099
4100    theta_ = 1.0e50;
4101    // now flip flop between spare arrays until reasonable theta
4102    tentativeTheta = CoinMax(10.0 * upperTheta, 1.0e-7);
4103
4104    // loops increasing tentative theta until can't go through
4105
4106    while (tentativeTheta < 1.0e22) {
4107      double thruThis = 0.0;
4108
4109      double bestPivot = acceptablePivot;
4110      int bestSequence = -1;
4111
4112      numberPossiblySwapped = numberColumns_;
4113      numberRemaining = 0;
4114
4115      upperTheta = 1.0e50;
4116
4117      spare = array[iFlip];
4118      index = indices[iFlip];
4119      spare2 = array[1 - iFlip];
4120      index2 = indices[1 - iFlip];
4121
4122      // try 3 different ways
4123      // 1 bias increase by ones with slightly wrong djs
4124      // 2 bias by all
4125      // 3 bias by all - tolerance
4126#define TRYBIAS 3
4127
4128      double increaseInThis = 0.0; //objective increase in this loop
4129
4130      for (i = 0; i < interesting[iFlip]; i++) {
4131        int iSequence = index[i];
4132        double alpha = spare[i];
4133        double oldValue = dj_[iSequence];
4134        double value = oldValue - tentativeTheta * alpha;
4135
4136        if (alpha < 0.0) {
4137          //at upper bound
4138          if (value > newTolerance) {
4139            double range = upper_[iSequence] - lower_[iSequence];
4140            thruThis -= range * alpha;
4141#if TRYBIAS == 1
4142            if (oldValue > 0.0)
4143              increaseInThis -= oldValue * range;
4144#elif TRYBIAS == 2
4145            increaseInThis -= oldValue * range;
4146#else
4147            increaseInThis -= (oldValue + dualTolerance_) * range;
4148#endif
4149            // goes on swapped list (also means candidates if too many)
4150            spare2[--numberPossiblySwapped] = alpha;
4151            index2[numberPossiblySwapped] = iSequence;
4152            if (fabs(alpha) > bestPivot) {
4153              bestPivot = fabs(alpha);
4154              bestSequence = numberPossiblySwapped;
4155            }
4156          } else {
4157            value = oldValue - upperTheta * alpha;
4158            if (value > newTolerance && -alpha >= acceptablePivot)
4159              upperTheta = (oldValue - newTolerance) / alpha;
4160            spare2[numberRemaining] = alpha;
4161            index2[numberRemaining++] = iSequence;
4162          }
4163        } else {
4164          // at lower bound
4165          if (value < -newTolerance) {
4166            double range = upper_[iSequence] - lower_[iSequence];
4167            thruThis += range * alpha;
4168            //?? is this correct - and should we look at good ones
4169#if TRYBIAS == 1
4170            if (oldValue < 0.0)
4171              increaseInThis += oldValue * range;
4172#elif TRYBIAS == 2
4173            increaseInThis += oldValue * range;
4174#else
4175            increaseInThis += (oldValue - dualTolerance_) * range;
4176#endif
4177            // goes on swapped list (also means candidates if too many)
4178            spare2[--numberPossiblySwapped] = alpha;
4179            index2[numberPossiblySwapped] = iSequence;
4180            if (fabs(alpha) > bestPivot) {
4181              bestPivot = fabs(alpha);
4182              bestSequence = numberPossiblySwapped;
4183            }
4184          } else {
4185            value = oldValue - upperTheta * alpha;
4186            if (value < -newTolerance && alpha >= acceptablePivot)
4187              upperTheta = (oldValue + newTolerance) / alpha;
4188            spare2[numberRemaining] = alpha;
4189            index2[numberRemaining++] = iSequence;
4190          }
4191        }
4192      }
4193      swapped[1 - iFlip] = numberPossiblySwapped;
4194      interesting[1 - iFlip] = numberRemaining;
4195      marker[1 - iFlip][0] = CoinMax(marker[1 - iFlip][0], numberRemaining);
4196      marker[1 - iFlip][1] = CoinMin(marker[1 - iFlip][1], numberPossiblySwapped);
4197
4198      double check = fabs(totalThru + thruThis);
4200      check += 1.0e-8 + 1.0e-10 * check;
4201      if (check >= fabs(dualOut_) || increaseInObjective + increaseInThis < 0.0) {
4202        // We should be pivoting in this batch
4203        // so compress down to this lot
4204        numberRemaining = 0;
4205        for (i = numberColumns_ - 1; i >= swapped[1 - iFlip]; i--) {
4206          spare[numberRemaining] = spare2[i];
4207          index[numberRemaining++] = index2[i];
4208        }
4209        interesting[iFlip] = numberRemaining;
4210        int iTry;
4211#define MAXTRY 100
4212        // first get ratio with tolerance
4213        for (iTry = 0; iTry < MAXTRY; iTry++) {
4214
4215          upperTheta = 1.0e50;
4216          numberPossiblySwapped = numberColumns_;
4217          numberRemaining = 0;
4218
4219          increaseInThis = 0.0; //objective increase in this loop
4220
4221          thruThis = 0.0;
4222
4223          spare = array[iFlip];
4224          index = indices[iFlip];
4225          spare2 = array[1 - iFlip];
4226          index2 = indices[1 - iFlip];
4227          for (i = 0; i < interesting[iFlip]; i++) {
4228            int iSequence = index[i];
4229            double alpha = spare[i];
4230            double oldValue = dj_[iSequence];
4231            double value = oldValue - upperTheta * alpha;
4232
4233            if (alpha < 0.0) {
4234              //at upper bound
4235              if (value > newTolerance) {
4236                if (-alpha >= acceptablePivot) {
4237                  upperTheta = (oldValue - newTolerance) / alpha;
4238                }
4239              }
4240            } else {
4241              // at lower bound
4242              if (value < -newTolerance) {
4243                if (alpha >= acceptablePivot) {
4244                  upperTheta = (oldValue + newTolerance) / alpha;
4245                }
4246              }
4247            }
4248          }
4249          bestPivot = acceptablePivot;
4250          sequenceIn_ = -1;
4251#ifdef DUBIOUS_WEIGHTS
4252          double bestWeight = COIN_DBL_MAX;
4253#endif
4254          double largestPivot = acceptablePivot;
4255          // now choose largest and sum all ones which will go through
4256          //printf("XX it %d number %d\n",numberIterations_,interesting[iFlip]);
4257          // Sum of bad small pivots
4258#ifdef MORE_CAREFUL
4261#endif
4262          // Make sure upperTheta will work (-O2 and above gives problems)
4263          upperTheta *= 1.0000000001;
4264          for (i = 0; i < interesting[iFlip]; i++) {
4265            int iSequence = index[i];
4266            double alpha = spare[i];
4267            double value = dj_[iSequence] - upperTheta * alpha;
4269
4271
4272            if (alpha < 0.0) {
4273              //at upper bound
4274              if (value >= 0.0) {
4276#if TRYBIAS == 1
4278#elif TRYBIAS == 2
4280#else
4281                badDj = -dj_[iSequence] - dualTolerance_;
4282#endif
4283              }
4284            } else {
4285              // at lower bound
4286              if (value <= 0.0) {
4288#if TRYBIAS == 1
4290#elif TRYBIAS == 2
4292#else
4293                badDj = dj_[iSequence] - dualTolerance_;
4294#endif
4295              }
4296            }
4298              // add to list of remaining
4299              spare2[numberRemaining] = alpha;
4300              index2[numberRemaining++] = iSequence;
4301            } else {
4302              // add to list of swapped
4303              spare2[--numberPossiblySwapped] = alpha;
4304              index2[numberPossiblySwapped] = iSequence;
4305              // select if largest pivot
4306              bool take = false;
4307              double absAlpha = fabs(alpha);
4308#ifdef DUBIOUS_WEIGHTS
4309              // User could do anything to break ties here
4310              double weight;
4311              if (dubiousWeights)
4312                weight = dubiousWeights[iSequence];
4313              else
4314                weight = 1.0;
4315              weight += randomNumberGenerator_.randomDouble() * 1.0e-2;
4316              if (absAlpha > 2.0 * bestPivot) {
4317                take = true;
4318              } else if (absAlpha > largestPivot) {
4319                // could multiply absAlpha and weight
4320                if (weight * bestPivot < bestWeight * absAlpha)
4321                  take = true;
4322              }
4323#else
4324              if (absAlpha > bestPivot)
4325                take = true;
4326#endif
4327#ifdef MORE_CAREFUL
4328              if (absAlpha < acceptablePivot && upperTheta < 1.0e20) {
4329                if (alpha < 0.0) {
4330                  //at upper bound
4331                  if (value > dualTolerance_) {
4332                    double gap = upper_[iSequence] - lower_[iSequence];
4333                    if (gap < 1.0e20)
4334                      sumBadPivots += value * gap;
4335                    else
4337                    //printf("bad %d alpha %g dj at upper %g\n",
4338                    //     iSequence,alpha,value);
4339                  }
4340                } else {
4341                  //at lower bound
4342                  if (value < -dualTolerance_) {
4343                    double gap = upper_[iSequence] - lower_[iSequence];
4344                    if (gap < 1.0e20)
4345                      sumBadPivots -= value * gap;
4346                    else
4348                    //printf("bad %d alpha %g dj at lower %g\n",
4349                    //     iSequence,alpha,value);
4350                  }
4351                }
4352              }
4353#endif
4354#ifdef FORCE_FOLLOW
4355              if (iSequence == force_in) {
4356                printf("taking %d - alpha %g best %g\n", force_in, absAlpha, largestPivot);
4357                take = true;
4358              }
4359#endif
4360              if (take) {
4361                sequenceIn_ = numberPossiblySwapped;
4362                bestPivot = absAlpha;
4363                theta_ = dj_[iSequence] / alpha;
4364                largestPivot = CoinMax(largestPivot, 0.5 * bestPivot);
4365#ifdef DUBIOUS_WEIGHTS
4366                bestWeight = weight;
4367#endif
4368                //printf(" taken seq %d alpha %g weight %d\n",
4369                //   iSequence,absAlpha,dubiousWeights[iSequence]);
4370              } else {
4371                //printf(" not taken seq %d alpha %g weight %d\n",
4372                //   iSequence,absAlpha,dubiousWeights[iSequence]);
4373              }
4374              double range = upper_[iSequence] - lower_[iSequence];
4375              thruThis += range * fabs(alpha);
4376              increaseInThis += badDj * range;
4377            }
4378          }
4379          marker[1 - iFlip][0] = CoinMax(marker[1 - iFlip][0], numberRemaining);
4380          marker[1 - iFlip][1] = CoinMin(marker[1 - iFlip][1], numberPossiblySwapped);
4381#ifdef MORE_CAREFUL
4382          // If we have done pivots and things look bad set alpha_ 0.0 to force factorization
4383          if (sumBadPivots > 1.0e4) {
4384            if (handler_->logLevel() > 1)
4385              *handler_ << "maybe forcing re-factorization - sum " << sumBadPivots << " " << factorization_->pivots() << " pivots" << CoinMessageEol;
4386            if (factorization_->pivots() > 3) {
4388              break;
4389            }
4390          }
4391#endif
4392          swapped[1 - iFlip] = numberPossiblySwapped;
4393          interesting[1 - iFlip] = numberRemaining;
4394          // If we stop now this will be increase in objective (I think)
4395          double increase = (fabs(dualOut_) - totalThru) * theta_;
4396          increase += increaseInObjective;
4397          if (theta_ < 0.0)
4398            thruThis += fabs(dualOut_); // force using this one
4399          if (increaseInObjective < 0.0 && increase < 0.0 && lastSequence >= 0) {
4400            // back
4401            // We may need to be more careful - we could do by
4402            // switch so we always do fine grained?
4403            bestPivot = 0.0;
4404          } else {
4406            totalThru += thruThis;
4407            increaseInObjective += increaseInThis;
4408          }
4409          if (bestPivot < 0.1 * bestEverPivot && bestEverPivot > 1.0e-6 && (bestPivot < 1.0e-3 || totalThru * 2.0 > fabs(dualOut_))) {
4410            // back to previous one
4411            sequenceIn_ = lastSequence;
4412            // swap regions
4413            iFlip = 1 - iFlip;
4414            break;
4415          } else if (sequenceIn_ == -1 && upperTheta > largeValue_) {
4416            if (lastPivot > acceptablePivot) {
4417              // back to previous one
4418              sequenceIn_ = lastSequence;
4419              // swap regions
4420              iFlip = 1 - iFlip;
4421            } else {
4422              // can only get here if all pivots too small
4423            }
4424            break;
4425          } else if (totalThru >= fabs(dualOut_)) {
4426            modifyCosts = true; // fine grain - we can modify costs
4427            break; // no point trying another loop
4428          } else {
4429            lastSequence = sequenceIn_;
4430            if (bestPivot > bestEverPivot)
4431              bestEverPivot = bestPivot;
4432            iFlip = 1 - iFlip;
4433            modifyCosts = true; // fine grain - we can modify costs
4434          }
4435        }
4436        if (iTry == MAXTRY)
4437          iFlip = 1 - iFlip; // flip back
4438        break;
4439      } else {
4440        // skip this lot
4441        if (bestPivot > 1.0e-3 || bestPivot > bestEverPivot) {
4442          bestEverPivot = bestPivot;
4443          lastSequence = bestSequence;
4444        } else {
4445          // keep old swapped
4446          CoinMemcpyN(array[iFlip] + swapped[iFlip],
4447            numberColumns_ - swapped[iFlip], array[1 - iFlip] + swapped[iFlip]);
4448          CoinMemcpyN(indices[iFlip] + swapped[iFlip],
4449            numberColumns_ - swapped[iFlip], indices[1 - iFlip] + swapped[iFlip]);
4450          marker[1 - iFlip][1] = CoinMin(marker[1 - iFlip][1], swapped[iFlip]);
4451          swapped[1 - iFlip] = swapped[iFlip];
4452        }
4453        increaseInObjective += increaseInThis;
4454        iFlip = 1 - iFlip; // swap regions
4455        tentativeTheta = 2.0 * upperTheta;
4456        totalThru += thruThis;
4457      }
4458    }
4459
4460    // can get here without sequenceIn_ set but with lastSequence
4461    if (sequenceIn_ < 0 && lastSequence >= 0) {
4462      // back to previous one
4463      sequenceIn_ = lastSequence;
4464      // swap regions
4465      iFlip = 1 - iFlip;
4466    }
4467
4468#define MINIMUMTHETA 1.0e-18
4469    // Movement should be minimum for anti-degeneracy - unless
4470    // fixed variable out
4471    double minimumTheta;
4472    if (upperOut_ > lowerOut_)
4473      minimumTheta = MINIMUMTHETA;
4474    else
4475      minimumTheta = 0.0;
4476    if (sequenceIn_ >= 0) {
4477      // at this stage sequenceIn_ is just pointer into index array
4478      // flip just so we can use iFlip
4479      iFlip = 1 - iFlip;
4480      spare = array[iFlip];
4481      index = indices[iFlip];
4482      double oldValue;
4483      alpha_ = spare[sequenceIn_];
4484      sequenceIn_ = indices[iFlip][sequenceIn_];
4485      oldValue = dj_[sequenceIn_];
4486      theta_ = CoinMax(oldValue / alpha_, 0.0);
4487      if (theta_ < minimumTheta && fabs(alpha_) < 1.0e5 && 1) {
4488        // can't pivot to zero
4489#if 0
4490                    if (oldValue - minimumTheta*alpha_ >= -dualTolerance_) {
4491                         theta_ = minimumTheta;
4492                    } else if (oldValue - minimumTheta*alpha_ >= -newTolerance) {
4493                         theta_ = minimumTheta;
4494                         thisIncrease = true;
4495                    } else {
4496                         theta_ = CoinMax((oldValue + newTolerance) / alpha_, 0.0);
4497                         thisIncrease = true;
4498                    }
4499#else
4500        theta_ = minimumTheta;
4501#endif
4502      }
4503      // may need to adjust costs so all dual feasible AND pivoted is exactly 0
4504      //int costOffset = numberRows_+numberColumns_;
4505      if (modifyCosts && !badSumPivots) {
4506        int i;
4507        for (i = numberColumns_ - 1; i >= swapped[iFlip]; i--) {
4508          int iSequence = index[i];
4509          double alpha = spare[i];
4510          double value = dj_[iSequence] - theta_ * alpha;
4511
4512          // can't be free here
4513
4514          if (alpha < 0.0) {
4515            //at upper bound
4516            if (value > dualTolerance_) {
4517              //thisIncrease = true;
4518#if CLP_CAN_HAVE_ZERO_OBJ < 2
4519#define MODIFYCOST 2
4520#endif
4521#if MODIFYCOST
4522              // modify cost to hit new tolerance
4523              double modification = alpha * theta_ - dj_[iSequence]
4524                + newTolerance;
4525              if ((specialOptions_ & (2048 + 4096 + 16384)) != 0) {
4526                if ((specialOptions_ & 16384) != 0) {
4527                  if (fabs(modification) < 1.0e-8)
4528                    modification = 0.0;
4529                } else if ((specialOptions_ & 2048) != 0) {
4530                  if (fabs(modification) < 1.0e-10)
4531                    modification = 0.0;
4532                } else {
4533                  if (fabs(modification) < 1.0e-12)
4534                    modification = 0.0;
4535                }
4536              }
4537              dj_[iSequence] += modification;
4538              cost_[iSequence] += modification;
4539              if (modification)
4540                numberChanged_++; // Say changed costs
4541                //cost_[iSequence+costOffset] += modification; // save change
4542#endif
4543            }
4544          } else {
4545            // at lower bound
4546            if (-value > dualTolerance_) {
4547              //thisIncrease = true;
4548#if MODIFYCOST
4549              // modify cost to hit new tolerance
4550              double modification = alpha * theta_ - dj_[iSequence]
4551                - newTolerance;
4552              //modification = CoinMax(modification,-dualTolerance_);
4553              //assert (fabs(modification)<1.0e-7);
4554              if ((specialOptions_ & (2048 + 4096)) != 0) {
4555                if ((specialOptions_ & 2048) != 0) {
4556                  if (fabs(modification) < 1.0e-10)
4557                    modification = 0.0;
4558                } else {
4559                  if (fabs(modification) < 1.0e-12)
4560                    modification = 0.0;
4561                }
4562              }
4563              dj_[iSequence] += modification;
4564              cost_[iSequence] += modification;
4565              if (modification)
4566                numberChanged_++; // Say changed costs
4567                //cost_[iSequence+costOffset] += modification; // save change
4568#endif
4569            }
4570          }
4571        }
4572      }
4573    }
4574  }
4575
4576#ifdef MORE_CAREFUL
4577  // If we have done pivots and things look bad set alpha_ 0.0 to force factorization
4578  if ((badSumPivots || fabs(theta_ * badFree) > 10.0 * dualTolerance_) && factorization_->pivots()) {
4579    if (handler_->logLevel() > 1)
4580      *handler_ << "forcing re-factorization" << CoinMessageEol;
4582    sequenceIn_ = -1;
4583    acceptablePivot_ = -acceptablePivot_;
4584  }
4585#endif
4586  if (sequenceIn_ >= 0) {
4587    lowerIn_ = lower_[sequenceIn_];
4588    upperIn_ = upper_[sequenceIn_];
4589    valueIn_ = solution_[sequenceIn_];
4590    dualIn_ = dj_[sequenceIn_];
4591
4592    if (numberTimesOptimal_) {
4593      // can we adjust cost back closer to original
4595    }
4596#if MODIFYCOST > 1
4597    // modify cost to hit zero exactly
4598    // so (dualIn_+modification)==theta_*alpha_
4599    double modification = theta_ * alpha_ - dualIn_;
4600    // But should not move objective too much ??
4601#define DONT_MOVE_OBJECTIVE
4602#ifdef DONT_MOVE_OBJECTIVE
4603    double moveObjective = fabs(modification * solution_[sequenceIn_]);
4604    double smallMove = CoinMax(fabs(objectiveValue_), 1.0e-3);
4605    if (moveObjective > smallMove) {
4606      if (handler_->logLevel() > 1)
4607        printf("would move objective by %g - original mod %g sol value %g\n", moveObjective,
4608          modification, solution_[sequenceIn_]);
4609      modification *= smallMove / moveObjective;
4610    }
4611#endif
4613      modification = 0.0;
4614    if ((specialOptions_ & (2048 + 4096)) != 0) {
4615      if ((specialOptions_ & 16384) != 0) {
4616        // in fast dual
4617        if (fabs(modification) < 1.0e-7)
4618          modification = 0.0;
4619      } else if ((specialOptions_ & 2048) != 0) {
4620        if (fabs(modification) < 1.0e-10)
4621          modification = 0.0;
4622      } else {
4623        if (fabs(modification) < 1.0e-12)
4624          modification = 0.0;
4625      }
4626    }
4627    dualIn_ += modification;
4628    dj_[sequenceIn_] = dualIn_;
4629    cost_[sequenceIn_] += modification;
4630    if (modification)
4631      numberChanged_++; // Say changed costs
4632      //int costOffset = numberRows_+numberColumns_;
4633      //cost_[sequenceIn_+costOffset] += modification; // save change
4634      //assert (fabs(modification)<1.0e-6);
4635#ifdef CLP_DEBUG
4636    if ((handler_->logLevel() & 32) && fabs(modification) > 1.0e-15)
4637      printf("exact %d new cost %g, change %g\n", sequenceIn_,
4638        cost_[sequenceIn_], modification);
4639#endif
4640#endif
4641
4642    if (alpha_ < 0.0) {
4643      // as if from upper bound
4644      directionIn_ = -1;
4645      upperIn_ = valueIn_;
4646    } else {
4647      // as if from lower bound
4648      directionIn_ = 1;
4649      lowerIn_ = valueIn_;
4650    }
4651    if (fabs(alpha_) < 1.0e-6) {
4652      // need bestPossible
4653      const double *work;
4654      int number;
4655      const int *which;
4656      const double *reducedCost;
4657      double tentativeTheta = 1.0e15;
4658      //double upperTheta = 1.0e31;
4659      bestPossible = 0.0;
4660      //double multiplier[] = { -1.0, 1.0 };
4661      double dualT = -dualTolerance_;
4662      int nSections = 2;
4664      for (int iSection = 0; iSection < nSections; iSection++) {
4665        if (!iSection) {
4666          work = rowArray->denseVector();
4667          number = rowArray->getNumElements();
4668          which = rowArray->getIndices();
4669          reducedCost = rowReducedCost_;
4671        } else {
4672          work = columnArray->denseVector();
4673          number = columnArray->getNumElements();
4674          which = columnArray->getIndices();
4675          reducedCost = reducedCostWork_;
4677        }
4678        for (i = 0; i < number; i++) {
4679          int iSequence = which[i];
4680          double alpha;
4681          double oldValue;
4682          double value;
4683          double mult = 1.0;
4684          switch (getStatus(iSequence + addSequence)) {
4685
4686          case basic:
4687          case ClpSimplex::isFixed:
4688            break;
4689          case isFree:
4690          case superBasic:
4691            alpha = work[i];
4692            bestPossible = CoinMax(bestPossible, fabs(alpha));
4693            break;
4694          case atUpperBound:
4695            mult = -1.0;
4696          case atLowerBound:
4697            alpha = work[i] * mult;
4698            if (alpha > 0.0) {
4699              oldValue = reducedCost[iSequence] * mult;
4700              value = oldValue - tentativeTheta * alpha;
4701              if (value < dualT) {
4702                bestPossible = CoinMax(bestPossible, alpha);
4703              }
4704            }
4705            break;
4706          }
4707        }
4708      }
4709    } else {
4710      bestPossible = fabs(alpha_);
4711    }
4712  } else {
4713    // no pivot
4714    bestPossible = 0.0;
4715    alpha_ = 0.0;
4716  }
4717  //if (thisIncrease)
4718  //dualTolerance_+= 1.0e-6*dblParam_[ClpDualTolerance];
4719
4720  // clear arrays
4721
4722  for (i = 0; i < 2; i++) {
4723    CoinZeroN(array[i], marker[i][0]);
4724    CoinZeroN(array[i] + marker[i][1], numberColumns_ - marker[i][1]);
4725  }
4726  return bestPossible;
4727}
4728#ifdef CLP_ALL_ONE_FILE
4729#undef MAXTRY
4730#endif
4731/* Checks if tentative optimal actually means unbounded
4732   Returns -3 if not, 2 if is unbounded */
4733int ClpSimplexDual::checkUnbounded(CoinIndexedVector *ray,
4734  CoinIndexedVector *spare,
4735  double changeCost)
4736{
4737  int status = 2; // say unbounded
4738  factorization_->updateColumn(spare, ray);
4739  // get reduced cost
4740  int i;
4741  int number = ray->getNumElements();
4742  int *index = ray->getIndices();
4743  double *array = ray->denseVector();
4744  for (i = 0; i < number; i++) {
4745    int iRow = index[i];
4746    int iPivot = pivotVariable_[iRow];
4747    changeCost -= cost(iPivot) * array[iRow];
4748  }
4749  double way;
4750  if (changeCost > 0.0) {
4751    //try going down
4752    way = 1.0;
4753  } else if (changeCost < 0.0) {
4754    //try going up
4755    way = -1.0;
4756  } else {
4757#ifdef CLP_DEBUG
4758    printf("can't decide on up or down\n");
4759#endif
4760    way = 0.0;
4761    status = -3;
4762  }
4763  double movement = 1.0e10 * way; // some largish number
4764  double zeroTolerance = 1.0e-14 * dualBound_;
4765  for (i = 0; i < number; i++) {
4766    int iRow = index[i];
4767    int iPivot = pivotVariable_[iRow];
4768    double arrayValue = array[iRow];
4769    if (fabs(arrayValue) < zeroTolerance)
4770      arrayValue = 0.0;
4771    double newValue = solution(iPivot) + movement * arrayValue;
4772    if (newValue > upper(iPivot) + primalTolerance_ || newValue < lower(iPivot) - primalTolerance_)
4773      status = -3; // not unbounded
4774  }
4775  if (status == 2) {
4776    // create ray
4777    delete[] ray_;
4778    ray_ = new double[numberColumns_];
4779    CoinZeroN(ray_, numberColumns_);
4780    for (i = 0; i < number; i++) {
4781      int iRow = index[i];
4782      int iPivot = pivotVariable_[iRow];
4783      double arrayValue = array[iRow];
4784      if (iPivot < numberColumns_ && fabs(arrayValue) >= zeroTolerance)
4785        ray_[iPivot] = way * array[iRow];
4786    }
4787  }
4788  ray->clear();
4789  return status;
4790}
4791//static int count_status=0;
4792//static double obj_status=0.0;
4793//static int check_status=123456789;//41754;
4794//static int count_alpha=0;
4795/* Checks if finished.  Updates status */
4796void ClpSimplexDual::statusOfProblemInDual(int &lastCleaned, int type,
4797  double *givenDuals, ClpDataSave &saveData,
4798  int ifValuesPass)
4799{
4800#ifdef CLP_INVESTIGATE_SERIAL
4801  if (z_thinks > 0 && z_thinks < 2)
4802    z_thinks += 2;
4803#endif
4804  bool arraysNotCreated = (type == 0);
4805  // If lots of iterations then adjust costs if large ones
4806  if (numberIterations_ > 4 * (numberRows_ + numberColumns_) && objectiveScale_ == 1.0) {
4807    double largest = 0.0;
4808    for (int i = 0; i < numberRows_; i++) {
4809      int iColumn = pivotVariable_[i];
4810      largest = CoinMax(largest, fabs(cost_[iColumn]));
4811    }
4812    if (largest > 1.0e6) {
4813      objectiveScale_ = 1.0e6 / largest;
4814      for (int i = 0; i < numberRows_ + numberColumns_; i++)
4815        cost_[i] *= objectiveScale_;
4816    }
4817  }
4818  int numberPivots = factorization_->pivots();
4819  double realDualInfeasibilities = 0.0;
4820  if (type == 2) {
4821    if (alphaAccuracy_ != -1.0)
4822      alphaAccuracy_ = -2.0;
4823    // trouble - restore solution
4824    CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4825    CoinMemcpyN(savedSolution_ + numberColumns_,
4826      numberRows_, rowActivityWork_);
4827    CoinMemcpyN(savedSolution_,
4828      numberColumns_, columnActivityWork_);
4829    // restore extra stuff
4830    int dummy;
4831    matrix_->generalExpanded(this, 6, dummy);
4832    forceFactorization_ = 1; // a bit drastic but ..
4833    changeMade_++; // say something changed
4834    // get correct bounds on all variables
4835    resetFakeBounds(0);
4836  }
4837  int tentativeStatus = problemStatus_;
4838  double changeCost;
4839  bool unflagVariables = true;
4840  bool weightsSaved = false;
4841  bool weightsSaved2 = numberIterations_ && !numberPrimalInfeasibilities_;
4842  int dontFactorizePivots = dontFactorizePivots_;
4843  if (type == 3) {
4844    type = 1;
4845    dontFactorizePivots = 1;
4846  }
4847  if (alphaAccuracy_ < 0.0 || !numberPivots || alphaAccuracy_ > 1.0e4 || numberPivots > 20) {
4848    if (problemStatus_ > -3 || numberPivots > dontFactorizePivots) {
4849      // factorize
4850      // later on we will need to recover from singularities
4851      // also we could skip if first time
4852      // save dual weights
4853      dualRowPivot_->saveWeights(this, 1);
4854      weightsSaved = true;
4855      if (type) {
4856        // is factorization okay?
4857        if (internalFactorize(1)) {
4858          // no - restore previous basis
4859          unflagVariables = false;
4860          assert(type == 1);
4861          changeMade_++; // say something changed
4862          // Keep any flagged variables
4863          int i;
4864          for (i = 0; i < numberRows_ + numberColumns_; i++) {
4865            if (flagged(i))
4866              saveStatus_[i] |= 64; //say flagged
4867          }
4868          CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4869          CoinMemcpyN(savedSolution_ + numberColumns_,
4870            numberRows_, rowActivityWork_);
4871          CoinMemcpyN(savedSolution_,
4872            numberColumns_, columnActivityWork_);
4873          // restore extra stuff
4874          int dummy;
4875          matrix_->generalExpanded(this, 6, dummy);
4876          // get correct bounds on all variables
4877          resetFakeBounds(1);
4878          // need to reject something
4879          char x = isColumn(sequenceOut_) ? 'C' : 'R';
4880          handler_->message(CLP_SIMPLEX_FLAG, messages_)
4881            << x << sequenceWithin(sequenceOut_)
4882            << CoinMessageEol;
4883#ifdef COIN_DEVELOP
4884          printf("flag d\n");
4885#endif
4886          setFlagged(sequenceOut_);
4888
4889          // Go to safe
4890          // not here - as can make more singular
4891          //factorization_->pivotTolerance(0.99);
4892          forceFactorization_ = 1; // a bit drastic but ..
4893          type = 2;
4894          //assert (internalFactorize(1)==0);
4895          if (internalFactorize(1)) {
4896            CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4897            CoinMemcpyN(savedSolution_ + numberColumns_,
4898              numberRows_, rowActivityWork_);
4899            CoinMemcpyN(savedSolution_,
4900              numberColumns_, columnActivityWork_);
4901            // restore extra stuff
4902            int dummy;
4903            matrix_->generalExpanded(this, 6, dummy);
4904            // debug
4905            int returnCode = internalFactorize(1);
4906            while (returnCode) {
4907              // ouch
4908              // switch off dense
4909              int saveDense = factorization_->denseThreshold();
4910              factorization_->setDenseThreshold(0);
4911              // Go to safe
4912              factorization_->pivotTolerance(0.99);
4913              // make sure will do safe factorization
4914              pivotVariable_[0] = -1;
4915              returnCode = internalFactorize(2);
4916              factorization_->setDenseThreshold(saveDense);
4917            }
4918            // get correct bounds on all variables
4919            resetFakeBounds(1);
4920          }
4921        }
4922      }
4923      if (problemStatus_ != -4 || numberPivots > 10)
4924        problemStatus_ = -3;
4925    }
4926  } else {
4927    //printf("testing with accuracy of %g and status of %d\n",alphaAccuracy_,problemStatus_);
4928    //count_alpha++;
4929    //if ((count_alpha%5000)==0)
4930    //printf("count alpha %d\n",count_alpha);
4931  }
4932  if (progress_.infeasibility_[0] < 1.0e-1 && primalTolerance_ == 1.0e-7 && progress_.iterationNumber_[0] > 0 && progress_.iterationNumber_[CLP_PROGRESS - 1] - progress_.iterationNumber_[0] > 25) {
4933    // default - so user did not set
4934    int iP;
4935    double minAverage = COIN_DBL_MAX;
4936    double maxAverage = 0.0;
4937    for (iP = 0; iP < CLP_PROGRESS; iP++) {
4938      int n = progress_.numberInfeasibilities_[iP];
4939      if (!n) {
4940        break;
4941      } else {
4942        double average = progress_.infeasibility_[iP];
4943        if (average > 0.1)
4944          break;
4945        average /= static_cast< double >(n);
4946        minAverage = CoinMin(minAverage, average);
4947        maxAverage = CoinMax(maxAverage, average);
4948      }
4949    }
4950    if (iP == CLP_PROGRESS && minAverage < 1.0e-5 && maxAverage < 1.0e-3) {
4951      // change tolerance
4952#if CBC_USEFUL_PRINTING > 0
4953      printf("CCchanging tolerance\n");
4954#endif
4955      primalTolerance_ = 1.0e-6;
4956      minimumPrimalTolerance_ = primalTolerance_;
4957      dblParam_[ClpPrimalTolerance] = 1.0e-6;
4958      moreSpecialOptions_ |= 4194304;
4959    }
4960  }
4961  // at this stage status is -3 or -4 if looks infeasible
4962  // get primal and dual solutions
4963#if 0
4964     {
4965          int numberTotal = numberRows_ + numberColumns_;
4966          double * saveSol = CoinCopyOfArray(solution_, numberTotal);
4967          double * saveDj = CoinCopyOfArray(dj_, numberTotal);
4968          double tolerance = type ? 1.0e-4 : 1.0e-8;
4969          // always if values pass
4970          double saveObj = objectiveValue_;
4971          double sumPrimal = sumPrimalInfeasibilities_;
4972          int numberPrimal = numberPrimalInfeasibilities_;
4973          double sumDual = sumDualInfeasibilities_;
4974          int numberDual = numberDualInfeasibilities_;
4975          gutsOfSolution(givenDuals, NULL);
4976          int j;
4977          double largestPrimal = tolerance;
4978          int iPrimal = -1;
4979          for (j = 0; j < numberTotal; j++) {
4980               double difference = solution_[j] - saveSol[j];
4981               if (fabs(difference) > largestPrimal) {
4982                    iPrimal = j;
4983                    largestPrimal = fabs(difference);
4984               }
4985          }
4986          double largestDual = tolerance;
4987          int iDual = -1;
4988          for (j = 0; j < numberTotal; j++) {
4989               double difference = dj_[j] - saveDj[j];
4990               if (fabs(difference) > largestDual && upper_[j] > lower_[j]) {
4991                    iDual = j;
4992                    largestDual = fabs(difference);
4993               }
4994          }
4995          if (!type) {
4996               if (fabs(saveObj - objectiveValue_) > 1.0e-5 ||
4997                         numberPrimal != numberPrimalInfeasibilities_ || numberPrimal != 1 ||
4998                         fabs(sumPrimal - sumPrimalInfeasibilities_) > 1.0e-5 || iPrimal >= 0 ||
4999                         numberDual != numberDualInfeasibilities_ || numberDual != 0 ||
5000                         fabs(sumDual - sumDualInfeasibilities_) > 1.0e-5 || iDual >= 0)
5001                    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",
5002                           type, numberIterations_, numberPivots,
5003                           numberPrimal, numberPrimalInfeasibilities_, sumPrimal, sumPrimalInfeasibilities_,
5004                           largestPrimal, iPrimal,
5005                           numberDual, numberDualInfeasibilities_, sumDual, sumDualInfeasibilities_,
5006                           largestDual, iDual,
5007                           saveObj, objectiveValue_);
5008          } else {
5009               if (fabs(saveObj - objectiveValue_) > 1.0e-5 ||
5010                         numberPrimalInfeasibilities_ || iPrimal >= 0 ||
5011                         numberDualInfeasibilities_ || iDual >= 0)
5012                    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",
5013                           type, numberIterations_, numberPivots,
5014                           numberPrimal, numberPrimalInfeasibilities_, sumPrimal, sumPrimalInfeasibilities_,
5015                           largestPrimal, iPrimal,
5016                           numberDual, numberDualInfeasibilities_, sumDual, sumDualInfeasibilities_,
5017                           largestDual, iDual,
5018                           saveObj, objectiveValue_);
5019          }
5020          delete [] saveSol;
5021          delete [] saveDj;
5022     }
5023#else
5024  if (type || ifValuesPass)
5025    gutsOfSolution(givenDuals, NULL);
5026#endif
5027  // If bad accuracy treat as singular
5028  if ((largestPrimalError_ > 1.0e15 || largestDualError_ > 1.0e15) && numberIterations_) {
5029    // restore previous basis
5030    unflagVariables = false;
5031    changeMade_++; // say something changed
5032    // Keep any flagged variables
5033    int i;
5034    for (i = 0; i < numberRows_ + numberColumns_; i++) {
5035      if (flagged(i))
5036        saveStatus_[i] |= 64; //say flagged
5037    }
5038    CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
5039    CoinMemcpyN(savedSolution_ + numberColumns_,
5040      numberRows_, rowActivityWork_);
5041    CoinMemcpyN(savedSolution_,
5042      numberColumns_, columnActivityWork_);
5043    // restore extra stuff
5044    int dummy;
5045    matrix_->generalExpanded(this, 6, dummy);
5046    // get correct bounds on all variables
5047    resetFakeBounds(1);
5048    // need to reject something
5049    int rejectedVariable = sequenceOut_;
5050    if (flagged(rejectedVariable)) {
5051      rejectedVariable = -1;
5052      for (int i = 0; i < numberRows_; i++) {
5053        int iSequence = pivotVariable_[i];
5054        if (!flagged(iSequence)) {
5055          rejectedVariable = iSequence;
5056          break;
5057        }
5058      }
5059      if (rejectedVariable < 0) {
5060        if (handler_->logLevel() > 1)
5061          printf("real trouble at line %d of ClpSimplexDual\n",
5062            __LINE__);
5063        problemStatus_ = 10;
5064        return;
5065      }
5066    }
5067    char x = isColumn(rejectedVariable) ? 'C' : 'R';
5068    handler_->message(CLP_SIMPLEX_FLAG, messages_)
5069      << x << sequenceWithin(rejectedVariable)
5070      << CoinMessageEol;
5071#ifdef COIN_DEVELOP
5072    printf("flag e\n");
5073#endif
5074    setFlagged(rejectedVariable);
5076
5077    // Go to safer
5078    double newTolerance = CoinMin(1.1 * factorization_->pivotTolerance(), 0.99);
5079    factorization_->pivotTolerance(newTolerance);
5080    forceFactorization_ = 1; // a bit drastic but ..
5081    if (alphaAccuracy_ != -1.0)
5082      alphaAccuracy_ = -2.0;
5083    type = 2;
5084    //assert (internalFactorize(1)==0);
5085    if (internalFactorize(1)) {
5086      CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
5087      CoinMemcpyN(savedSolution_ + numberColumns_,
5088        numberRows_, rowActivityWork_);
5089      CoinMemcpyN(savedSolution_,
5090        numberColumns_, columnActivityWork_);
5091      // restore extra stuff
5092      int dummy;
5093      matrix_->generalExpanded(this, 6, dummy);
5094      // debug
5095      int returnCode = internalFactorize(1);
5096      while (returnCode) {
5097        // ouch
5098        // switch off dense
5099        int saveDense = factorization_->denseThreshold();
5100        factorization_->setDenseThreshold(0);
5101        // Go to safe
5102        factorization_->pivotTolerance(0.99);
5103        // make sure will do safe factorization
5104        pivotVariable_[0] = -1;
5105        returnCode = internalFactorize(2);
5106        factorization_->setDenseThreshold(saveDense);
5107      }
5108      // get correct bounds on all variables
5109      resetFakeBounds(1);
5110    }
5111    // get primal and dual solutions
5112    gutsOfSolution(givenDuals, NULL);
5113  } else if (goodAccuracy()) {
5114    // Can reduce tolerance
5115    double newTolerance = CoinMax(0.995 * factorization_->pivotTolerance(), saveData.pivotTolerance_);
5116    factorization_->pivotTolerance(newTolerance);
5117  }
5118  bestObjectiveValue_ = CoinMax(bestObjectiveValue_,
5119    objectiveValue_ - bestPossibleImprovement_);
5121  // Double check infeasibility if no action
5122  if (progress_.lastIterationNumber(0) == numberIterations_) {
5123    if (dualRowPivot_->looksOptimal()) {
5124      numberPrimalInfeasibilities_ = 0;
5125      sumPrimalInfeasibilities_ = 0.0;
5126    }
5127#if 1
5128  } else {
5129    double thisObj = objectiveValue_ - bestPossibleImprovement_;
5130#ifdef CLP_INVESTIGATE
5131    assert(bestPossibleImprovement_ > -1000.0 && objectiveValue_ > -1.0e100);
5132    if (bestPossibleImprovement_)
5133      printf("obj %g add in %g -> %g\n", objectiveValue_, bestPossibleImprovement_,
5134        thisObj);
5135#endif
5136    double lastObj = progress_.lastObjective(0);
5137#ifndef NDEBUG
5138#ifdef COIN_DEVELOP
5139    resetFakeBounds(-1);
5140#endif
5141#endif
5142#ifdef CLP_REPORT_PROGRESS
5143    ixxxxxx++;
5144    if (ixxxxxx >= ixxyyyy - 4 && ixxxxxx <= ixxyyyy) {
5145      char temp[20];
5146      sprintf(temp, "sol%d.out", ixxxxxx);
5147      printf("sol%d.out\n", ixxxxxx);
5148      FILE *fp = fopen(temp, "w");
5149      int nTotal = numberRows_ + numberColumns_;
5150      for (int i = 0; i < nTotal; i++)
5151        fprintf(fp, "%d %d %g %g %g %g %g\n",
5152          i, status_[i], lower_[i], solution_[i], upper_[i], cost_[i], dj_[i]);
5153      fclose(fp);
5154    }
5155#endif
5156    if (!ifValuesPass && firstFree_ < 0) {
5157      double testTol = 5.0e-3;
5158      if (progress_.timesFlagged() > 10) {
5159        testTol *= pow(2.0, progress_.timesFlagged() - 8);
5160      } else if (progress_.timesFlagged() > 5) {
5161        testTol *= 5.0;
5162      }
5163      if (lastObj > thisObj + testTol * (fabs(thisObj) + fabs(lastObj)) + testTol) {
5164        int maxFactor = factorization_->maximumPivots();
5165        if ((specialOptions_ & 1048576) == 0) {
5166          if (progress_.timesFlagged() > 10)
5168          if (maxFactor > 10 - 9) {
5169#ifdef COIN_DEVELOP
5170            printf("lastobj %g thisobj %g\n", lastObj, thisObj);
5171#endif
5172            //if (forceFactorization_<0)
5173            //forceFactorization_= maxFactor;
5174            //forceFactorization_ = CoinMax(1,(forceFactorization_>>1));
5175            if ((progressFlag_ & 4) == 0 && lastObj < thisObj + 1.0e4 && largestPrimalError_ < 1.0e2) {
5176              // Just save costs
5177              // save extra copy of cost_
5178              int nTotal = numberRows_ + numberColumns_;
5179              double *temp = new double[2 * nTotal];
5180              memcpy(temp, cost_, nTotal * sizeof(double));
5181              memcpy(temp + nTotal, cost_, nTotal * sizeof(double));
5182              delete[] cost_;
5183              cost_ = temp;
5184              objectiveWork_ = cost_;
5185              rowObjectiveWork_ = cost_ + numberColumns_;
5186              progressFlag_ |= 4;
5187            } else {
5188              forceFactorization_ = 1;
5189#ifdef COIN_DEVELOP
5190              printf("Reducing factorization frequency - bad backwards\n");
5191#endif
5192#if 1
5193              unflagVariables = false;
5194              changeMade_++; // say something changed
5195              int nTotal = numberRows_ + numberColumns_;
5196              CoinMemcpyN(saveStatus_, nTotal, status_);
5197              CoinMemcpyN(savedSolution_ + numberColumns_,
5198                numberRows_, rowActivityWork_);
5199              CoinMemcpyN(savedSolution_,
5200                numberColumns_, columnActivityWork_);
5201              if ((progressFlag_ & 4) == 0) {
5202                // save extra copy of cost_
5203                double *temp = new double[2 * nTotal];
5204                memcpy(temp, cost_, nTotal * sizeof(double));
5205                memcpy(temp + nTotal, cost_, nTotal * sizeof(double));
5206                delete[] cost_;
5207                cost_ = temp;
5208                objectiveWork_ = cost_;
5209                rowObjectiveWork_ = cost_ + numberColumns_;
5210                progressFlag_ |= 4;
5211              } else {
5212                memcpy(cost_, cost_ + nTotal, nTotal * sizeof(double));
5213              }
5214              // restore extra stuff
5215              int dummy;
5216              matrix_->generalExpanded(this, 6, dummy);
5217              double pivotTolerance = factorization_->pivotTolerance();
5218              if (pivotTolerance < 0.2)
5219                factorization_->pivotTolerance(0.2);
5220              else if (progress_.timesFlagged() > 2)
5221                factorization_->pivotTolerance(CoinMin(pivotTolerance * 1.1, 0.99));
5222              if (alphaAccuracy_ != -1.0)
5223                alphaAccuracy_ = -2.0;
5224              if (internalFactorize(1)) {
5225                CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
5226                CoinMemcpyN(savedSolution_ + numberColumns_,
5227                  numberRows_, rowActivityWork_);
5228                CoinMemcpyN(savedSolution_,
5229                  numberColumns_, columnActivityWork_);
5230                // restore extra stuff
5231                int dummy;
5232                matrix_->generalExpanded(this, 6, dummy);
5233                // debug
5234                int returnCode = internalFactorize(1);
5235                while (returnCode) {
5236                  // ouch
5237                  // switch off dense
5238                  int saveDense = factorization_->denseThreshold();
5239                  factorization_->setDenseThreshold(0);
5240                  // Go to safe
5241                  factorization_->pivotTolerance(0.99);
5242                  // make sure will do safe factorization
5243                  pivotVariable_[0] = -1;
5244                  returnCode = internalFactorize(2);
5245                  factorization_->setDenseThreshold(saveDense);
5246                }
5247              }
5248              resetFakeBounds(0);
5249              type = 2; // so will restore weights
5250              // get primal and dual solutions
5251              gutsOfSolution(givenDuals, NULL);
5252              if (numberPivots < 2) {
5253                // need to reject something
5254                char x = isColumn(sequenceOut_) ? 'C' : 'R';
5255                handler_->message(CLP_SIMPLEX_FLAG, messages_)
5256                  << x << sequenceWithin(sequenceOut_)
5257                  << CoinMessageEol;
5258#ifdef COIN_DEVELOP
5259                printf("flag d\n");
5260#endif
5261                setFlagged(sequenceOut_);
5263                progress_.incrementTimesFlagged();
5264              }
5265              if (numberPivots < 10)
5267#ifdef COIN_DEVELOP
5268              printf("obj now %g\n", objectiveValue_);
5269#endif
5270              progress_.modifyObjective(objectiveValue_
5271                - bestPossibleImprovement_);
5272#endif
5273            }
5274          }
5275        } else {
5276          // in fast dual give up
5277#ifdef COIN_DEVELOP
5278          printf("In fast dual?\n");
5279#endif
5280          problemStatus_ = 3;
5281        }
5282      } else if (lastObj < thisObj - 1.0e-5 * CoinMax(fabs(thisObj), fabs(lastObj)) - 1.0e-3) {
5283        numberTimesOptimal_ = 0;
5284      }
5285    }
5286#endif
5287  }
5288  // Up tolerance if looks a bit odd
5289  if (numberIterations_ > CoinMax(1000, numberRows_ >> 4) && (specialOptions_ & 64) != 0) {
5290    if (sumPrimalInfeasibilities_ && sumPrimalInfeasibilities_ < 1.0e5) {
5291      int backIteration = progress_.lastIterationNumber(CLP_PROGRESS - 1);
5292      if (backIteration > 0 && numberIterations_ - backIteration < 9 * CLP_PROGRESS) {
5293        if (factorization_->pivotTolerance() < 0.9) {
5294          // up tolerance
5295          factorization_->pivotTolerance(CoinMin(factorization_->pivotTolerance() * 1.05 + 0.02, 0.91));
5296          //printf("tol now %g\n",factorization_->pivotTolerance());
5297          progress_.clearIterationNumbers();
5298        }
5299      }
5300    }
5301  }
5302  // Check if looping
5303  int loop;
5304  if (!givenDuals && type != 2)
5305    loop = progress_.looping();
5306  else
5307    loop = -1;
5308  if (progress_.reallyBadTimes() > 10) {
5309    problemStatus_ = 10; // instead - try other algorithm
5310#if COIN_DEVELOP > 2
5311    printf("returning at %d\n", __LINE__);
5312#endif
5313  }
5314  int situationChanged = 0;
5315  if (loop >= 0) {
5316    problemStatus_ = loop; //exit if in loop
5317    if (!problemStatus_) {
5318      // declaring victory
5319      numberPrimalInfeasibilities_ = 0;
5320      sumPrimalInfeasibilities_ = 0.0;
5321    } else {
5322      problemStatus_ = 10; // instead - try other algorithm
5323#if COIN_DEVELOP > 2
5324      printf("returning at %d\n", __LINE__);
5325#endif
5326    }
5327    return;
5328  } else if (loop < -1) {
5329    // something may have changed
5330    gutsOfSolution(NULL, NULL);
5331    situationChanged = 1;
5332  }
5334  if ((progressFlag_ & 2) != 0) {
5335    situationChanged = 2;
5336  }
5337  progressFlag_ &= (~3); //reset progress flag
5338  if ((progressFlag_ & 4) != 0) {
5339    // save copy of cost_
5340    int nTotal = numberRows_ + numberColumns_;
5341    memcpy(cost_ + nTotal, cost_, nTotal * sizeof(double));
5342  }
5343  /*if (!numberIterations_&&sumDualInfeasibilities_)
5344       printf("OBJ %g sumPinf %g sumDinf %g\n",
5345        objectiveValue(),sumPrimalInfeasibilities_,
5346        sumDualInfeasibilities_);*/
5347  // mark as having gone optimal if looks like it
5348  if (!numberPrimalInfeasibilities_ && !numberDualInfeasibilities_)
5349    progressFlag_ |= 8;
5350  if (handler_->detail(CLP_SIMPLEX_STATUS, messages_) < 100) {
5351    handler_->message(CLP_SIMPLEX_STATUS, messages_)
5352      << numberIterations_ << objectiveValue();
5353    handler_->printing(sumPrimalInfeasibilities_ > 0.0)
5354      << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
5355    handler_->printing(sumDualInfeasibilities_ > 0.0)
5356      << sumDualInfeasibilities_ << numberDualInfeasibilities_;
5357    handler_->printing(numberDualInfeasibilitiesWithoutFree_
5358      < numberDualInfeasibilities_)
5359      << numberDualInfeasibilitiesWithoutFree_;
5360    handler_->message() << CoinMessageEol;
5361  }
5362#if 0
5363     count_status++;
5364     if (!numberIterations_)
5365       obj_status=-1.0e30;
5366     if (objectiveValue()<obj_status-0.01) {
5367       printf("Backward obj at %d from %g to %g\n",
5368              count_status,obj_status,objectiveValue());
5369     }
5370     obj_status=objectiveValue();
5371     if (count_status>=check_status-1) {
5372       printf("Trouble ahead - count_status %d\n",count_status);
5373     }
5374#endif
5375#if 0
5376     printf("IT %d %g %g(%d) %g(%d)\n",
5377            numberIterations_, objectiveValue(),
5378            sumPrimalInfeasibilities_, numberPrimalInfeasibilities_,
5379            sumDualInfeasibilities_, numberDualInfeasibilities_);
5380#endif
5381  double approximateObjective = objectiveValue_;
5382#ifdef CLP_REPORT_PROGRESS
5383  if (ixxxxxx >= ixxyyyy - 4 && ixxxxxx <= ixxyyyy) {
5384    char temp[20];
5385    sprintf(temp, "x_sol%d.out", ixxxxxx);
5386    FILE *fp = fopen(temp, "w");
5387    int nTotal = numberRows_ + numberColumns_;
5388    for (int i = 0; i < nTotal; i++)
5389      fprintf(fp, "%d %d %g %g %g %g %g\n",
5390        i, status_[i], lower_[i], solution_[i], upper_[i], cost_[i], dj_[i]);
5391    fclose(fp);
5392    if (ixxxxxx == ixxyyyy)
5393      exit(6);
5394  }
5395#endif
5396  realDualInfeasibilities = sumDualInfeasibilities_;
5397  double saveTolerance = dualTolerance_;
5398  // If we need to carry on cleaning variables
5399  if (!numberPrimalInfeasibilities_ && (specialOptions_ & 1024) != 0 && CLEAN_FIXED) {
5400    for (int iRow = 0; iRow < numberRows_; iRow++) {
5401      int iPivot = pivotVariable_[iRow];
5402      if (!flagged(iPivot) && pivoted(iPivot)) {
5403        // carry on
5404        numberPrimalInfeasibilities_ = -1;
5405        sumOfRelaxedPrimalInfeasibilities_ = 1.0;
5406        sumPrimalInfeasibilities_ = 1.0;
5407        break;
5408      }
5409    }
5410  }
5411  /* If we are primal feasible and any dual infeasibilities are on
5412        free variables then it is better to go to primal */
5413  if (!numberPrimalInfeasibilities_ && ((!numberDualInfeasibilitiesWithoutFree_ && numberDualInfeasibilities_) || (moreSpecialOptions_ & 16777216) != 0))
5414    problemStatus_ = 10;
5415  // dual bound coming in
5416  double saveDualBound = dualBound_;
5417  bool needCleanFake = false;
5418  while (problemStatus_ <= -3 && saveDualBound == dualBound_) {
5419    int cleanDuals = 0;
5420    if (situationChanged != 0)
5421      cleanDuals = 1;
5422    int numberChangedBounds = 0;
5423    int doOriginalTolerance = 0;
5424    if (lastCleaned == numberIterations_)
5425      doOriginalTolerance = 1;
5426    // check optimal
5427    // give code benefit of doubt
5428    if (sumOfRelaxedDualInfeasibilities_ == 0.0 && sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
5429      // say optimal (with these bounds etc)
5430      numberDualInfeasibilities_ = 0;
5431      sumDualInfeasibilities_ = 0.0;
5432      numberPrimalInfeasibilities_ = 0;
5433      sumPrimalInfeasibilities_ = 0.0;
5434    }
5435    //if (dualFeasible()||problemStatus_==-4||(primalFeasible()&&!numberDualInfeasibilitiesWithoutFree_)) {
5436    if (dualFeasible() || problemStatus_ == -4) {
5437      progress_.modifyObjective(objectiveValue_
5438        - bestPossibleImprovement_);
5439#ifdef COIN_DEVELOP
5440      if (sumDualInfeasibilities_ || bestPossibleImprovement_)
5441        printf("improve %g dualinf %g -> %g\n",
5442          bestPossibleImprovement_, sumDualInfeasibilities_,
5443          sumDualInfeasibilities_ * dualBound_);
5444#endif
5445      // see if cutoff reached
5446      double limit = 0.0;
5447      getDblParam(ClpDualObjectiveLimit, limit);
5448#if 0
5449               if(fabs(limit) < 1.0e30 && objectiveValue()*optimizationDirection_ >
5450                         limit + 1.0e-7 + 1.0e-8 * fabs(limit) && !numberAtFakeBound()) {
5451                    //looks infeasible on objective
5452                    if (perturbation_ == 101) {
5453                         cleanDuals = 1;
5454                         // Save costs
5455                         int numberTotal = numberRows_ + numberColumns_;
5456                         double * saveCost = CoinCopyOfArray(cost_, numberTotal);
5457                         // make sure fake bounds are back
5458                         changeBounds(1, NULL, changeCost);
5459                         createRim4(false);
5460                         // make sure duals are current
5461                         computeDuals(givenDuals);
5462                         checkDualSolution();
5463                         if(objectiveValue()*optimizationDirection_ >
5464                                   limit + 1.0e-7 + 1.0e-8 * fabs(limit) && !numberDualInfeasibilities_) {
5465                              perturbation_ = 102; // stop any perturbations
5466                              printf("cutoff test succeeded\n");
5467                         } else {
5468                              printf("cutoff test failed\n");
5469                              // put back
5470                              memcpy(cost_, saveCost, numberTotal * sizeof(double));
5471                              // make sure duals are current
5472                              computeDuals(givenDuals);
5473                              checkDualSolution();
5474                              progress_.modifyObjective(-COIN_DBL_MAX);
5475                              problemStatus_ = -1;
5476                         }
5477                         delete [] saveCost;
5478                    }
5479               }
5480#endif
5481      if (primalFeasible() && !givenDuals) {
5482        // may be optimal - or may be bounds are wrong
5483        handler_->message(CLP_DUAL_BOUNDS, messages_)
5484          << dualBound_
5485          << CoinMessageEol;
5486        // save solution in case unbounded
5487        double *saveColumnSolution = NULL;
5488        double *saveRowSolution = NULL;
5489        bool inCbc = (specialOptions_ & (0x01000000 | 16384)) != 0;
5490        if (!inCbc) {
5491          saveColumnSolution = CoinCopyOfArray(columnActivityWork_, numberColumns_);
5492          saveRowSolution = CoinCopyOfArray(rowActivityWork_, numberRows_);
5493        }
5494#ifndef COIN_MAX_DUAL_BOUND
5495#define COIN_MAX_DUAL_BOUND 1.0e20
5496#endif
5497        numberChangedBounds = (dualBound_ < COIN_MAX_DUAL_BOUND) ? changeBounds(0, rowArray_[3], changeCost) : 0;
5498        if (numberChangedBounds <= 0 && !numberDualInfeasibilities_) {
5499          //looks optimal - do we need to reset tolerance
5500          if (perturbation_ == 101) {
5501            perturbation_ = 102; // stop any perturbations
5502            cleanDuals = 1;
5503            // make sure fake bounds are back
5504            //computeObjectiveValue();
5505            changeBounds(1, NULL, changeCost);
5506            //computeObjectiveValue();
5507            createRim4(false);
5508            // make sure duals are current
5509            computeDuals(givenDuals);
5510            checkDualSolution();
5511            progress_.modifyObjective(-COIN_DBL_MAX);
5512#define DUAL_TRY_FASTER
5513#ifdef DUAL_TRY_FASTER
5514            if (numberDualInfeasibilities_) {
5515#endif
5516              numberChanged_ = 1; // force something to happen
5517              lastCleaned = numberIterations_ - 1;
5518#ifdef DUAL_TRY_FASTER
5519            } else {
5520              //double value = objectiveValue_;
5521              computeObjectiveValue(true);
5522              //printf("old %g new %g\n",value,objectiveValue_);
5523              //numberChanged_=1;
5524            }
5525#endif
5526          }
5527          if (lastCleaned < numberIterations_ && numberTimesOptimal_ < 4 && (numberChanged_ || (specialOptions_ & 4096) == 0)) {
5528#if CLP_CAN_HAVE_ZERO_OBJ
5529            if ((specialOptions_ & 16777216) == 0) {
5530#endif
5531              doOriginalTolerance = 2;
5532              numberTimesOptimal_++;
5533              changeMade_++; // say something changed
5534              if (numberTimesOptimal_ == 1) {
5535                dualTolerance_ = dblParam_[ClpDualTolerance];
5536              } else {
5537                if (numberTimesOptimal_ == 2) {
5538                  // better to have small tolerance even if slower
5539                  factorization_->zeroTolerance(CoinMin(factorization_->zeroTolerance(), 1.0e-15));
5540                }
5541                dualTolerance_ = dblParam_[ClpDualTolerance];
5542                dualTolerance_ *= pow(2.0, numberTimesOptimal_ - 1);
5543              }
5544              cleanDuals = 2; // If nothing changed optimal else primal
5545#if CLP_CAN_HAVE_ZERO_OBJ
5546            } else {
5547              // no cost - skip checks
5548              problemStatus_ = 0;
5549            }
5550#endif
5551          } else {
5552            problemStatus_ = 0; // optimal
5553            if (lastCleaned < numberIterations_ && numberChanged_) {
5554              handler_->message(CLP_SIMPLEX_GIVINGUP, messages_)
5555                << CoinMessageEol;
5556            }
5557          }
5558        } else {
5559          cleanDuals = 1;
5560          if (doOriginalTolerance == 1) {
5561            // check unbounded
5562            // find a variable with bad dj
5563            int iSequence;
5564            int iChosen = -1;
5565            if (!inCbc) {
5566              double largest = 100.0 * primalTolerance_;
5567              for (iSequence = 0; iSequence < numberRows_ + numberColumns_;
5568                   iSequence++) {
5569                double djValue = dj_[iSequence];
5570                double originalLo = originalLower(iSequence);
5571                double originalUp = originalUpper(iSequence);
5572                if (fabs(djValue) > fabs(largest)) {
5573                  if (getStatus(iSequence) != basic) {
5574                    if (djValue > 0 && originalLo < -1.0e20) {
5575                      if (djValue > fabs(largest)) {
5576                        largest = djValue;
5577                        iChosen = iSequence;
5578                      }
5579                    } else if (djValue < 0 && originalUp > 1.0e20) {
5580                      if (-djValue > fabs(largest)) {
5581                        largest = djValue;
5582                        iChosen = iSequence;
5583                      }
5584                    }
5585                  }
5586                }
5587              }
5588            }
5589            if (iChosen >= 0) {
5590              int iSave = sequenceIn_;
5591              sequenceIn_ = iChosen;
5592              unpack(rowArray_[1]);
5593              sequenceIn_ = iSave;
5594              // if dual infeasibilities then must be free vector so add in dual
5595              if (numberDualInfeasibilities_) {
5596                if (fabs(changeCost) > 1.0e-5)
5597                  COIN_DETAIL_PRINT(printf("Odd free/unbounded combo