# source:stable/1.16/Clp/src/ClpSimplexDual.cpp@2181

Last change on this file since 2181 was 2181, checked in by forrest, 4 years ago

out accurate duals sometime

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